Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/boost/math/special_functions/beta.hpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 //  (C) Copyright John Maddock 2006.
0002 //  (C) Copyright Matt Borland 2024.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 
0007 #ifndef BOOST_MATH_SPECIAL_BETA_HPP
0008 #define BOOST_MATH_SPECIAL_BETA_HPP
0009 
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013 
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/tools/type_traits.hpp>
0016 #include <boost/math/tools/assert.hpp>
0017 #include <boost/math/tools/precision.hpp>
0018 #include <boost/math/tools/numeric_limits.hpp>
0019 #include <boost/math/tools/tuple.hpp>
0020 #include <boost/math/tools/promotion.hpp>
0021 #include <boost/math/tools/cstdint.hpp>
0022 #include <boost/math/special_functions/gamma.hpp>
0023 #include <boost/math/special_functions/erf.hpp>
0024 #include <boost/math/special_functions/log1p.hpp>
0025 #include <boost/math/special_functions/expm1.hpp>
0026 #include <boost/math/special_functions/trunc.hpp>
0027 #include <boost/math/special_functions/lanczos.hpp>
0028 #include <boost/math/policies/policy.hpp>
0029 #include <boost/math/policies/error_handling.hpp>
0030 #include <boost/math/constants/constants.hpp>
0031 #include <boost/math/special_functions/math_fwd.hpp>
0032 #include <boost/math/special_functions/binomial.hpp>
0033 #include <boost/math/special_functions/factorials.hpp>
0034 #include <boost/math/tools/roots.hpp>
0035 
0036 namespace boost{ namespace math{
0037 
0038 namespace detail{
0039 
0040 //
0041 // Implementation of Beta(a,b) using the Lanczos approximation:
0042 //
0043 template <class T, class Lanczos, class Policy>
0044 BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
0045 {
0046    BOOST_MATH_STD_USING  // for ADL of std names
0047 
0048    if(a <= 0)
0049       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);
0050    if(b <= 0)
0051       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);
0052 
0053    T result;  // LCOV_EXCL_LINE
0054 
0055    T prefix = 1;
0056    T c = a + b;
0057 
0058    // Special cases:
0059    if((c == a) && (b < tools::epsilon<T>()))
0060       return 1 / b;
0061    else if((c == b) && (a < tools::epsilon<T>()))
0062       return 1 / a;
0063    if(b == 1)
0064       return 1/a;
0065    else if(a == 1)
0066       return 1/b;
0067    else if(c < tools::epsilon<T>())
0068    {
0069       result = c / a;
0070       result /= b;
0071       return result;
0072    }
0073 
0074    /*
0075    //
0076    // This code appears to be no longer necessary: it was
0077    // used to offset errors introduced from the Lanczos
0078    // approximation, but the current Lanczos approximations
0079    // are sufficiently accurate for all z that we can ditch
0080    // this.  It remains in the file for future reference...
0081    //
0082    // If a or b are less than 1, shift to greater than 1:
0083    if(a < 1)
0084    {
0085       prefix *= c / a;
0086       c += 1;
0087       a += 1;
0088    }
0089    if(b < 1)
0090    {
0091       prefix *= c / b;
0092       c += 1;
0093       b += 1;
0094    }
0095    */
0096 
0097    if(a < b)
0098    {
0099       BOOST_MATH_GPU_SAFE_SWAP(a, b);
0100    }
0101 
0102    // Lanczos calculation:
0103    T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
0104    T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
0105    T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
0106    result = Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c));
0107    T ambh = a - 0.5f - b;
0108    if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
0109    {
0110       // Special case where the base of the power term is close to 1
0111       // compute (1+x)^y instead:
0112       result *= exp(ambh * boost::math::log1p(-b / cgh, pol));
0113    }
0114    else
0115    {
0116       result *= pow(agh / cgh, a - T(0.5) - b);
0117    }
0118    if(cgh > 1e10f)
0119       // this avoids possible overflow, but appears to be marginally less accurate:
0120       result *= pow((agh / cgh) * (bgh / cgh), b);
0121    else
0122       result *= pow((agh * bgh) / (cgh * cgh), b);
0123    result *= sqrt(boost::math::constants::e<T>() / bgh);
0124 
0125    // If a and b were originally less than 1 we need to scale the result:
0126    result *= prefix;
0127 
0128    return result;
0129 } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&)
0130 
0131 //
0132 // Generic implementation of Beta(a,b) without Lanczos approximation support
0133 // (Caution this is slow!!!):
0134 //
0135 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0136 template <class T, class Policy>
0137 BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol)
0138 {
0139    BOOST_MATH_STD_USING
0140 
0141    if(a <= 0)
0142       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);
0143    if(b <= 0)
0144       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);
0145 
0146    const T c = a + b;
0147 
0148    // Special cases:
0149    if ((c == a) && (b < tools::epsilon<T>()))
0150       return 1 / b;
0151    else if ((c == b) && (a < tools::epsilon<T>()))
0152       return 1 / a;
0153    if (b == 1)
0154       return 1 / a;
0155    else if (a == 1)
0156       return 1 / b;
0157    else if (c < tools::epsilon<T>())
0158    {
0159       T result = c / a;
0160       result /= b;
0161       return result;
0162    }
0163 
0164    // Regular cases start here:
0165    const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
0166 
0167    long shift_a = 0;
0168    long shift_b = 0;
0169 
0170    if(a < min_sterling)
0171       shift_a = 1 + ltrunc(min_sterling - a);
0172    if(b < min_sterling)
0173       shift_b = 1 + ltrunc(min_sterling - b);
0174    long shift_c = shift_a + shift_b;
0175 
0176    if ((shift_a == 0) && (shift_b == 0))
0177    {
0178       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);
0179    }
0180    else if ((a < 1) && (b < 1))
0181    {
0182       return boost::math::tgamma(a, pol) * (boost::math::tgamma(b, pol) / boost::math::tgamma(c));
0183    }
0184    else if(a < 1)
0185       return boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol);
0186    else if(b < 1)
0187       return boost::math::tgamma(b, pol) * boost::math::tgamma_delta_ratio(a, b, pol);
0188    else
0189    {
0190       T result = beta_imp(T(a + shift_a), T(b + shift_b), l, pol);
0191       //
0192       // Recursion:
0193       //
0194       for (long i = 0; i < shift_c; ++i)
0195       {
0196          result *= c + i;
0197          if (i < shift_a)
0198             result /= a + i;
0199          if (i < shift_b)
0200             result /= b + i;
0201       }
0202       return result;
0203    }
0204 
0205 } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l)
0206 #endif
0207 
0208 //
0209 // Compute the leading power terms in the incomplete Beta:
0210 //
0211 // (x^a)(y^b)/Beta(a,b) when normalised, and
0212 // (x^a)(y^b) otherwise.
0213 //
0214 // Almost all of the error in the incomplete beta comes from this
0215 // function: particularly when a and b are large. Computing large
0216 // powers are *hard* though, and using logarithms just leads to
0217 // horrendous cancellation errors.
0218 //
0219 template <class T, class Lanczos, class Policy>
0220 BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a,
0221                         T b,
0222                         T x,
0223                         T y,
0224                         const Lanczos&,
0225                         bool normalised,
0226                         const Policy& pol,
0227                         T prefix = 1,
0228                         const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
0229 {
0230    BOOST_MATH_STD_USING
0231 
0232    if(!normalised)
0233    {
0234       // can we do better here?
0235       return pow(x, a) * pow(y, b);
0236    }
0237 
0238    T result;  // LCOV_EXCL_LINE
0239 
0240    T c = a + b;
0241 
0242    // combine power terms with Lanczos approximation:
0243    T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
0244    T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
0245    T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
0246    if ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
0247       result = 0;  // denominator overflows in this case
0248    else
0249       result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
0250    result *= prefix;
0251    // combine with the leftover terms from the Lanczos approximation:
0252    result *= sqrt(bgh / boost::math::constants::e<T>());
0253    result *= sqrt(agh / cgh);
0254 
0255    // l1 and l2 are the base of the exponents minus one:
0256    T l1 = (x * b - y * agh) / agh;
0257    T l2 = (y * a - x * bgh) / bgh;
0258    if((BOOST_MATH_GPU_SAFE_MIN(fabs(l1), fabs(l2)) < 0.2))
0259    {
0260       // when the base of the exponent is very near 1 we get really
0261       // gross errors unless extra care is taken:
0262       if((l1 * l2 > 0) || (BOOST_MATH_GPU_SAFE_MIN(a, b) < 1))
0263       {
0264          //
0265          // This first branch handles the simple cases where either:
0266          //
0267          // * The two power terms both go in the same direction
0268          // (towards zero or towards infinity).  In this case if either
0269          // term overflows or underflows, then the product of the two must
0270          // do so also.
0271          // *Alternatively if one exponent is less than one, then we
0272          // can't productively use it to eliminate overflow or underflow
0273          // from the other term.  Problems with spurious overflow/underflow
0274          // can't be ruled out in this case, but it is *very* unlikely
0275          // since one of the power terms will evaluate to a number close to 1.
0276          //
0277          if(fabs(l1) < 0.1)
0278          {
0279             result *= exp(a * boost::math::log1p(l1, pol));
0280             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0281          }
0282          else
0283          {
0284             result *= pow((x * cgh) / agh, a);
0285             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0286          }
0287          if(fabs(l2) < 0.1)
0288          {
0289             result *= exp(b * boost::math::log1p(l2, pol));
0290             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0291          }
0292          else
0293          {
0294             result *= pow((y * cgh) / bgh, b);
0295             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0296          }
0297       }
0298       else if(BOOST_MATH_GPU_SAFE_MAX(fabs(l1), fabs(l2)) < 0.5)
0299       {
0300          //
0301          // Both exponents are near one and both the exponents are
0302          // greater than one and further these two
0303          // power terms tend in opposite directions (one towards zero,
0304          // the other towards infinity), so we have to combine the terms
0305          // to avoid any risk of overflow or underflow.
0306          //
0307          // We do this by moving one power term inside the other, we have:
0308          //
0309          //    (1 + l1)^a * (1 + l2)^b
0310          //  = ((1 + l1)*(1 + l2)^(b/a))^a
0311          //  = (1 + l1 + l3 + l1*l3)^a   ;  l3 = (1 + l2)^(b/a) - 1
0312          //                                    = exp((b/a) * log(1 + l2)) - 1
0313          //
0314          // The tricky bit is deciding which term to move inside :-)
0315          // By preference we move the larger term inside, so that the
0316          // size of the largest exponent is reduced.  However, that can
0317          // only be done as long as l3 (see above) is also small.
0318          //
0319          bool small_a = a < b;
0320          T ratio = b / a;
0321          if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1)))
0322          {
0323             T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol);
0324             l3 = l1 + l3 + l3 * l1;
0325             l3 = a * boost::math::log1p(l3, pol);
0326             result *= exp(l3);
0327             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0328          }
0329          else
0330          {
0331             T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol);
0332             l3 = l2 + l3 + l3 * l2;
0333             l3 = b * boost::math::log1p(l3, pol);
0334             result *= exp(l3);
0335             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0336          }
0337       }
0338       else if(fabs(l1) < fabs(l2))
0339       {
0340          // First base near 1 only:
0341          T l = a * boost::math::log1p(l1, pol)
0342             + b * log((y * cgh) / bgh);
0343          if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
0344          {
0345             l += log(result);
0346             if(l >= tools::log_max_value<T>())
0347                return policies::raise_overflow_error<T>(function, nullptr, pol);  // LCOV_EXCL_LINE we can probably never get here, probably!
0348             result = exp(l);
0349          }
0350          else
0351             result *= exp(l);
0352          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0353       }
0354       else
0355       {
0356          // Second base near 1 only:
0357          T l = b * boost::math::log1p(l2, pol)
0358             + a * log((x * cgh) / agh);
0359          if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>()))
0360          {
0361             l += log(result);
0362             if(l >= tools::log_max_value<T>())
0363                return policies::raise_overflow_error<T>(function, nullptr, pol);  // LCOV_EXCL_LINE we can probably never get here, probably!
0364             result = exp(l);
0365          }
0366          else
0367             result *= exp(l);
0368          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0369       }
0370    }
0371    else
0372    {
0373       // general case:
0374       T b1 = (x * cgh) / agh;
0375       T b2 = (y * cgh) / bgh;
0376       l1 = a * log(b1);
0377       l2 = b * log(b2);
0378       BOOST_MATH_INSTRUMENT_VARIABLE(b1);
0379       BOOST_MATH_INSTRUMENT_VARIABLE(b2);
0380       BOOST_MATH_INSTRUMENT_VARIABLE(l1);
0381       BOOST_MATH_INSTRUMENT_VARIABLE(l2);
0382       if((l1 >= tools::log_max_value<T>())
0383          || (l1 <= tools::log_min_value<T>())
0384          || (l2 >= tools::log_max_value<T>())
0385          || (l2 <= tools::log_min_value<T>())
0386          )
0387       {
0388          // Oops, under/overflow, sidestep if we can:
0389          if(a < b)
0390          {
0391             T p1 = pow(b2, b / a);
0392             T l3 = (b1 != 0) && (p1 != 0) ? (a * (log(b1) + log(p1))) : tools::max_value<T>(); // arbitrary large value if the logs would fail!
0393             if((l3 < tools::log_max_value<T>())
0394                && (l3 > tools::log_min_value<T>()))
0395             {
0396                result *= pow(p1 * b1, a);
0397             }
0398             else
0399             {
0400                l2 += l1 + log(result);
0401                if(l2 >= tools::log_max_value<T>())
0402                   return policies::raise_overflow_error<T>(function, nullptr, pol);  // LCOV_EXCL_LINE we can probably never get here, probably!
0403                result = exp(l2);
0404             }
0405          }
0406          else
0407          {
0408             // This protects against spurious overflow in a/b:
0409             T p1 = (b1 < 1) && (b < 1) && (tools::max_value<T>() * b < a) ? static_cast<T>(0) : static_cast<T>(pow(b1, a / b));
0410             T l3 = (p1 != 0) && (b2 != 0) ? (log(p1) + log(b2)) * b : tools::max_value<T>();  // arbitrary large value if the logs would fail!
0411             if((l3 < tools::log_max_value<T>())
0412                && (l3 > tools::log_min_value<T>()))
0413             {
0414                result *= pow(p1 * b2, b);
0415             }
0416             else if(result != 0)  // we can elude the calculation below if we're already going to be zero
0417             {
0418                l2 += l1 + log(result);
0419                if(l2 >= tools::log_max_value<T>())
0420                   return policies::raise_overflow_error<T>(function, nullptr, pol);  // LCOV_EXCL_LINE we can probably never get here, probably!
0421                result = exp(l2);
0422             }
0423          }
0424          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0425       }
0426       else
0427       {
0428          // finally the normal case:
0429          result *= pow(b1, a) * pow(b2, b);
0430          BOOST_MATH_INSTRUMENT_VARIABLE(result);
0431       }
0432    }
0433 
0434    BOOST_MATH_INSTRUMENT_VARIABLE(result);
0435 
0436    if (0 == result)
0437    {
0438       if ((a > 1) && (x == 0))
0439          return result;  // true zero LCOV_EXCL_LINE we can probably never get here
0440       if ((b > 1) && (y == 0))
0441          return result; // true zero LCOV_EXCL_LINE we can probably never get here
0442       return boost::math::policies::raise_underflow_error<T>(function, nullptr, pol);
0443    }
0444 
0445    return result;
0446 }
0447 //
0448 // Compute the leading power terms in the incomplete Beta:
0449 //
0450 // (x^a)(y^b)/Beta(a,b) when normalised, and
0451 // (x^a)(y^b) otherwise.
0452 //
0453 // Almost all of the error in the incomplete beta comes from this
0454 // function: particularly when a and b are large. Computing large
0455 // powers are *hard* though, and using logarithms just leads to
0456 // horrendous cancellation errors.
0457 //
0458 // This version is generic, slow, and does not use the Lanczos approximation.
0459 //
0460 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0461 template <class T, class Policy>
0462 BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a,
0463                         T b,
0464                         T x,
0465                         T y,
0466                         const boost::math::lanczos::undefined_lanczos& l,
0467                         bool normalised,
0468                         const Policy& pol,
0469                         T prefix = 1,
0470                         const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)")
0471 {
0472    BOOST_MATH_STD_USING
0473 
0474    if(!normalised)
0475    {
0476       return prefix * pow(x, a) * pow(y, b);
0477    }
0478 
0479    T c = a + b;
0480 
0481    const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
0482 
0483    long shift_a = 0;
0484    long shift_b = 0;
0485 
0486    if (a < min_sterling)
0487       shift_a = 1 + ltrunc(min_sterling - a);
0488    if (b < min_sterling)
0489       shift_b = 1 + ltrunc(min_sterling - b);
0490 
0491    if ((shift_a == 0) && (shift_b == 0))
0492    {
0493       T power1, power2;
0494       bool need_logs = false;
0495       if (a < b)
0496       {
0497          BOOST_MATH_IF_CONSTEXPR(boost::math::numeric_limits<T>::has_infinity)
0498          {
0499             power1 = pow((x * y * c * c) / (a * b), a);
0500             power2 = pow((y * c) / b, b - a);
0501          }
0502          else
0503          {
0504             // We calculate these logs purely so we can check for overflow in the power functions
0505             T l1 = log((x * y * c * c) / (a * b));
0506             T l2 = log((y * c) / b);
0507             if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
0508             {
0509                power1 = pow((x * y * c * c) / (a * b), a);
0510                power2 = pow((y * c) / b, b - a);
0511             }
0512             else
0513             {
0514                need_logs = true;
0515             }
0516          }
0517       }
0518       else
0519       {
0520          BOOST_MATH_IF_CONSTEXPR(boost::math::numeric_limits<T>::has_infinity)
0521          {
0522             power1 = pow((x * y * c * c) / (a * b), b);
0523             power2 = pow((x * c) / a, a - b);
0524          }
0525          else
0526          {
0527             // We calculate these logs purely so we can check for overflow in the power functions
0528             T l1 = log((x * y * c * c) / (a * b)) * b;
0529             T l2 = log((x * c) / a) * (a - b);
0530             if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
0531             {
0532                power1 = pow((x * y * c * c) / (a * b), b);
0533                power2 = pow((x * c) / a, a - b);
0534             }
0535             else
0536                need_logs = true;
0537          }
0538       }
0539       BOOST_MATH_IF_CONSTEXPR(boost::math::numeric_limits<T>::has_infinity)
0540       {
0541          if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2))
0542          {
0543             need_logs = true;
0544          }
0545       }
0546       if (need_logs)
0547       {
0548          //
0549          // We want:
0550          //
0551          // (xc / a)^a (yc / b)^b
0552          //
0553          // But we know that one or other term will over / underflow and combining the logs will be next to useless as that will cause significant cancellation.
0554          // If we assume b > a and express z ^ b as(z ^ b / a) ^ a with z = (yc / b) then we can move one power term inside the other :
0555          //
0556          // ((xc / a) * (yc / b)^(b / a))^a
0557          //
0558          // However, we're not quite there yet, as the term being exponentiated is quite likely to be close to unity, so let:
0559          //
0560          // xc / a = 1 + (xb - ya) / a
0561          //
0562          // analogously let :
0563          //
0564          // 1 + p = (yc / b) ^ (b / a) = 1 + expm1((b / a) * log1p((ya - xb) / b))
0565          //
0566          // so putting the two together we have :
0567          //
0568          // exp(a * log1p((xb - ya) / a + p + p(xb - ya) / a))
0569          //
0570          // Analogously, when a > b we can just swap all the terms around.
0571          //
0572          // Finally, there are a few cases (x or y is unity) when the above logic can't be used
0573          // or where there is no logarithmic cancellation and accuracy is better just using
0574          // the regular formula:
0575          //
0576          T xc_a = x * c / a;
0577          T yc_b = y * c / b;
0578          if ((x == 1) || (y == 1) || (fabs(xc_a - 1) > 0.25) || (fabs(yc_b - 1) > 0.25))
0579          {
0580             // The above logic fails, the result is almost certainly zero:
0581             power1 = exp(log(xc_a) * a + log(yc_b) * b);
0582             power2 = 1;
0583          }
0584          else if (b > a)
0585          {
0586             T p = boost::math::expm1((b / a) * boost::math::log1p((y * a - x * b) / b));
0587             power1 = exp(a * boost::math::log1p((x * b - y * a) / a + p * (x * c / a)));
0588             power2 = 1;
0589          }
0590          else
0591          {
0592             T p = boost::math::expm1((a / b) * boost::math::log1p((x * b - y * a) / a));
0593             power1 = exp(b * boost::math::log1p((y * a - x * b) / b + p * (y * c / b)));
0594             power2 = 1;
0595          }
0596       }
0597       return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
0598    }
0599 
0600    T power1 = pow(x, a);
0601    T power2 = pow(y, b);
0602    T bet = beta_imp(a, b, l, pol);
0603 
0604    if(!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2) || !(boost::math::isnormal)(bet))
0605    {
0606       int shift_c = shift_a + shift_b;
0607       T result = ibeta_power_terms(T(a + shift_a), T(b + shift_b), x, y, l, normalised, pol, prefix);
0608       if ((boost::math::isnormal)(result))
0609       {
0610          for (int i = 0; i < shift_c; ++i)
0611          {
0612             result /= c + i;
0613                if (i < shift_a)
0614                {
0615                   result *= a + i;
0616                      result /= x;
0617                }
0618             if (i < shift_b)
0619             {
0620                result *= b + i;
0621                result /= y;
0622             }
0623          }
0624          return prefix * result;
0625       }
0626       else
0627       {
0628          T log_result = log(x) * a + log(y) * b + log(prefix);
0629          if ((boost::math::isnormal)(bet))
0630             log_result -= log(bet);
0631          else
0632             log_result += boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol);
0633          return exp(log_result);
0634       }
0635    }
0636    return prefix * power1 * (power2 / bet);
0637 }
0638 
0639 #endif
0640 //
0641 // Series approximation to the incomplete beta:
0642 //
0643 template <class T>
0644 struct ibeta_series_t
0645 {
0646    typedef T result_type;
0647    BOOST_MATH_GPU_ENABLED ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {}
0648    BOOST_MATH_GPU_ENABLED T operator()()
0649    {
0650       T r = result / apn;
0651       apn += 1;
0652       result *= poch * x / n;
0653       ++n;
0654       poch += 1;
0655       return r;
0656    }
0657 private:
0658    T result, x, apn, poch;
0659    int n;
0660 };
0661 
0662 template <class T, class Lanczos, class Policy>
0663 BOOST_MATH_GPU_ENABLED T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
0664 {
0665    BOOST_MATH_STD_USING
0666 
0667    T result;
0668 
0669    BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
0670 
0671    if(normalised)
0672    {
0673       T c = a + b;
0674 
0675       // incomplete beta power term, combined with the Lanczos approximation:
0676       T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
0677       T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
0678       T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
0679       if ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
0680          result = 0;  // denorms cause overflow in the Lanzos series, result will be zero anyway
0681       else
0682          result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
0683 
0684       if (!(boost::math::isfinite)(result))
0685          result = 0;  // LCOV_EXCL_LINE we can probably never get here, covered already above?
0686 
0687       T l1 = log(cgh / bgh) * (b - 0.5f);
0688       T l2 = log(x * cgh / agh) * a;
0689       //
0690       // Check for over/underflow in the power terms:
0691       //
0692       if((l1 > tools::log_min_value<T>())
0693          && (l1 < tools::log_max_value<T>())
0694          && (l2 > tools::log_min_value<T>())
0695          && (l2 < tools::log_max_value<T>()))
0696       {
0697          if(a * b < bgh * 10)
0698             result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol));
0699          else
0700             result *= pow(cgh / bgh, T(b - T(0.5)));
0701          result *= pow(x * cgh / agh, a);
0702          result *= sqrt(agh / boost::math::constants::e<T>());
0703 
0704          if(p_derivative)
0705          {
0706             *p_derivative = result * pow(y, b);
0707             BOOST_MATH_ASSERT(*p_derivative >= 0);
0708          }
0709       }
0710       else
0711       {
0712          //
0713          // Oh dear, we need logs, and this *will* cancel:
0714          //
0715          if (result != 0)  // elude calculation when result will be zero.
0716          {
0717             result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
0718             if (p_derivative)
0719                *p_derivative = exp(result + b * log(y));
0720             result = exp(result);
0721          }
0722       }
0723    }
0724    else
0725    {
0726       // Non-normalised, just compute the power:
0727       result = pow(x, a);
0728    }
0729    if(result < tools::min_value<T>())
0730       return s0; // Safeguard: series can't cope with denorms.
0731    ibeta_series_t<T> s(a, b, x, result);
0732    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0733    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
0734    policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol);
0735    return result;
0736 }
0737 //
0738 // Incomplete Beta series again, this time without Lanczos support:
0739 //
0740 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0741 template <class T, class Policy>
0742 BOOST_MATH_GPU_ENABLED 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)
0743 {
0744    BOOST_MATH_STD_USING
0745 
0746    T result;
0747    BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
0748 
0749    if(normalised)
0750    {
0751       const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
0752 
0753       long shift_a = 0;
0754       long shift_b = 0;
0755 
0756       if (a < min_sterling)
0757          shift_a = 1 + ltrunc(min_sterling - a);
0758       if (b < min_sterling)
0759          shift_b = 1 + ltrunc(min_sterling - b);
0760 
0761       T c = a + b;
0762 
0763       if ((shift_a == 0) && (shift_b == 0))
0764       {
0765          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));
0766       }
0767       else if ((a < 1) && (b > 1))
0768          result = pow(x, a) / (boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol));
0769       else
0770       {
0771          T power = pow(x, a);
0772          T bet = beta_imp(a, b, l, pol);
0773          if (!(boost::math::isnormal)(power) || !(boost::math::isnormal)(bet))
0774          {
0775             result = exp(a * log(x) + boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol));
0776          }
0777          else
0778             result = power / bet;
0779       }
0780       if(p_derivative)
0781       {
0782          *p_derivative = result * pow(y, b);
0783          BOOST_MATH_ASSERT(*p_derivative >= 0);
0784       }
0785    }
0786    else
0787    {
0788       // Non-normalised, just compute the power:
0789       result = pow(x, a);
0790    }
0791    if(result < tools::min_value<T>())
0792       return s0; // Safeguard: series can't cope with denorms.
0793    ibeta_series_t<T> s(a, b, x, result);
0794    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0795    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0);
0796    policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol);
0797    return result;
0798 }
0799 #endif
0800 //
0801 // Continued fraction for the incomplete beta:
0802 //
0803 template <class T>
0804 struct ibeta_fraction2_t
0805 {
0806    typedef boost::math::pair<T, T> result_type;
0807 
0808    BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {}
0809 
0810    BOOST_MATH_GPU_ENABLED result_type operator()()
0811    {
0812       T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x;
0813       T denom = (a + 2 * m - 1);
0814       aN /= denom * denom;
0815 
0816       T bN = static_cast<T>(m);
0817       bN += (m * (b - m) * x) / (a + 2*m - 1);
0818       bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1);
0819 
0820       ++m;
0821 
0822       return boost::math::make_pair(aN, bN);
0823    }
0824 
0825 private:
0826    T a, b, x, y;
0827    int m;
0828 };
0829 //
0830 // Evaluate the incomplete beta via the continued fraction representation:
0831 //
0832 template <class T, class Policy>
0833 BOOST_MATH_GPU_ENABLED inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
0834 {
0835    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
0836    BOOST_MATH_STD_USING
0837    T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
0838    if(p_derivative)
0839    {
0840       *p_derivative = result;
0841       BOOST_MATH_ASSERT(*p_derivative >= 0);
0842    }
0843    if(result == 0)
0844       return result;
0845 
0846    ibeta_fraction2_t<T> f(a, b, x, y);
0847    T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>());
0848    BOOST_MATH_INSTRUMENT_VARIABLE(fract);
0849    BOOST_MATH_INSTRUMENT_VARIABLE(result);
0850    return result / fract;
0851 }
0852 //
0853 // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x):
0854 //
0855 template <class T, class Policy>
0856 BOOST_MATH_GPU_ENABLED T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
0857 {
0858    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
0859 
0860    BOOST_MATH_INSTRUMENT_VARIABLE(k);
0861 
0862    T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
0863    if(p_derivative)
0864    {
0865       *p_derivative = prefix;
0866       BOOST_MATH_ASSERT(*p_derivative >= 0);
0867    }
0868    prefix /= a;
0869    if(prefix == 0)
0870       return prefix;
0871    T sum = 1;
0872    T term = 1;
0873    // series summation from 0 to k-1:
0874    for(int i = 0; i < k-1; ++i)
0875    {
0876       term *= (a+b+i) * x / (a+i+1);
0877       sum += term;
0878    }
0879    prefix *= sum;
0880 
0881    return prefix;
0882 }
0883 
0884 //
0885 // This function is only needed for the non-regular incomplete beta,
0886 // it computes the delta in:
0887 // beta(a,b,x) = prefix + delta * beta(a+k,b,x)
0888 // it is currently only called for small k.
0889 //
0890 template <class T>
0891 BOOST_MATH_GPU_ENABLED inline T rising_factorial_ratio(T a, T b, int k)
0892 {
0893    // calculate:
0894    // (a)(a+1)(a+2)...(a+k-1)
0895    // _______________________
0896    // (b)(b+1)(b+2)...(b+k-1)
0897 
0898    // This is only called with small k, for large k
0899    // it is grossly inefficient, do not use outside it's
0900    // intended purpose!!!
0901    BOOST_MATH_INSTRUMENT_VARIABLE(k);
0902    BOOST_MATH_ASSERT(k > 0);
0903 
0904    T result = 1;
0905    for(int i = 0; i < k; ++i)
0906       result *= (a+i) / (b+i);
0907    return result;
0908 }
0909 //
0910 // Routine for a > 15, b < 1
0911 //
0912 // Begin by figuring out how large our table of Pn's should be,
0913 // quoted accuracies are "guesstimates" based on empirical observation.
0914 // Note that the table size should never exceed the size of our
0915 // tables of factorials.
0916 //
0917 template <class T>
0918 struct Pn_size
0919 {
0920    // This is likely to be enough for ~35-50 digit accuracy
0921    // but it's hard to quantify exactly:
0922    #ifndef BOOST_MATH_HAS_NVRTC
0923    static constexpr unsigned value =
0924       ::boost::math::max_factorial<T>::value >= 100 ? 50
0925    : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<double>::value ? 30
0926    : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value ? 15 : 1;
0927    static_assert(::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value, "Type does not provide for 35-50 digits of accuracy.");
0928    #else
0929    static constexpr unsigned value = 0; // Will never be called
0930    #endif
0931 };
0932 template <>
0933 struct Pn_size<float>
0934 {
0935    static constexpr unsigned value = 15; // ~8-15 digit accuracy
0936 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0937    static_assert(::boost::math::max_factorial<float>::value >= 30, "Type does not provide for 8-15 digits of accuracy.");
0938 #endif
0939 };
0940 template <>
0941 struct Pn_size<double>
0942 {
0943    static constexpr unsigned value = 30; // 16-20 digit accuracy
0944 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0945    static_assert(::boost::math::max_factorial<double>::value >= 60, "Type does not provide for 16-20 digits of accuracy.");
0946 #endif
0947 };
0948 template <>
0949 struct Pn_size<long double>
0950 {
0951    static constexpr unsigned value = 50; // ~35-50 digit accuracy
0952 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0953    static_assert(::boost::math::max_factorial<long double>::value >= 100, "Type does not provide for ~35-50 digits of accuracy");
0954 #endif
0955 };
0956 
0957 template <class T, class Policy>
0958 BOOST_MATH_GPU_ENABLED T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised)
0959 {
0960    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
0961    BOOST_MATH_STD_USING
0962    //
0963    // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6.
0964    //
0965    // Some values we'll need later, these are Eq 9.1:
0966    //
0967    T bm1 = b - 1;
0968    T t = a + bm1 / 2;
0969    T lx, u;  // LCOV_EXCL_LINE
0970    if(y < 0.35)
0971       lx = boost::math::log1p(-y, pol);
0972    else
0973       lx = log(x);
0974    u = -t * lx;
0975    // and from from 9.2:
0976    T prefix;  // LCOV_EXCL_LINE
0977    T h = regularised_gamma_prefix(b, u, pol, lanczos_type());
0978    if(h <= tools::min_value<T>())
0979       return s0;
0980    if(normalised)
0981    {
0982       prefix = h / boost::math::tgamma_delta_ratio(a, b, pol);
0983       prefix /= pow(t, b);
0984    }
0985    else
0986    {
0987       prefix = full_igamma_prefix(b, u, pol) / pow(t, b);
0988    }
0989    prefix *= mult;
0990    //
0991    // now we need the quantity Pn, unfortunately this is computed
0992    // recursively, and requires a full history of all the previous values
0993    // so no choice but to declare a big table and hope it's big enough...
0994    //
0995    T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 };  // see 9.3.
0996    //
0997    // Now an initial value for J, see 9.6:
0998    //
0999    T j = boost::math::gamma_q(b, u, pol) / h;
1000    //
1001    // Now we can start to pull things together and evaluate the sum in Eq 9:
1002    //
1003    T sum = s0 + prefix * j;  // Value at N = 0
1004    // some variables we'll need:
1005    unsigned tnp1 = 1; // 2*N+1
1006    T lx2 = lx / 2;
1007    lx2 *= lx2;
1008    T lxp = 1;
1009    T t4 = 4 * t * t;
1010    T b2n = b;
1011 
1012    for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n)
1013    {
1014       /*
1015       // debugging code, enable this if you want to determine whether
1016       // the table of Pn's is large enough...
1017       //
1018       static int max_count = 2;
1019       if(n > max_count)
1020       {
1021          max_count = n;
1022          std::cerr << "Max iterations in BGRAT was " << n << std::endl;
1023       }
1024       */
1025       //
1026       // begin by evaluating the next Pn from Eq 9.4:
1027       //
1028       tnp1 += 2;
1029       p[n] = 0;
1030       T mbn = b - n;
1031       unsigned tmp1 = 3;
1032       for(unsigned m = 1; m < n; ++m)
1033       {
1034          mbn = m * b - n;
1035          p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1);
1036          tmp1 += 2;
1037       }
1038       p[n] /= n;
1039       p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1);
1040       //
1041       // Now we want Jn from Jn-1 using Eq 9.6:
1042       //
1043       j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4;
1044       lxp *= lx2;
1045       b2n += 2;
1046       //
1047       // pull it together with Eq 9:
1048       //
1049       T r = prefix * p[n] * j;
1050       sum += r;
1051       // r is always small:
1052       BOOST_MATH_ASSERT(tools::max_value<T>() * tools::epsilon<T>() > fabs(r));
1053       if(fabs(r / tools::epsilon<T>()) < fabs(sum))
1054          break;
1055    }
1056    return sum;
1057 } // 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)
1058 
1059 //
1060 // For integer arguments we can relate the incomplete beta to the
1061 // complement of the binomial distribution cdf and use this finite sum.
1062 //
1063 template <class T, class Policy>
1064 BOOST_MATH_GPU_ENABLED T binomial_ccdf(T n, T k, T x, T y, const Policy& pol)
1065 {
1066    BOOST_MATH_STD_USING // ADL of std names
1067 
1068    T result = pow(x, n);
1069 
1070    if(result > tools::min_value<T>())
1071    {
1072       T term = result;
1073       for(unsigned i = itrunc(T(n - 1)); i > k; --i)
1074       {
1075          term *= ((i + 1) * y) / ((n - i) * x);
1076          result += term;
1077       }
1078    }
1079    else
1080    {
1081       // First term underflows so we need to start at the mode of the
1082       // distribution and work outwards:
1083       int start = itrunc(n * x);
1084       if(start <= k + 1)
1085          start = itrunc(k + 2);
1086       result = static_cast<T>(pow(x, T(start)) * pow(y, n - T(start)) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start), pol));
1087       if(result == 0)
1088       {
1089          // OK, starting slightly above the mode didn't work,
1090          // we'll have to sum the terms the old fashioned way.
1091          // Very hard to get here, possibly only when exponent
1092          // range is very limited (as with type float):
1093          // LCOV_EXCL_START
1094          for(unsigned i = start - 1; i > k; --i)
1095          {
1096             result += static_cast<T>(pow(x, static_cast<T>(i)) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i), pol));
1097          }
1098          // LCOV_EXCL_STOP
1099       }
1100       else
1101       {
1102          T term = result;
1103          T start_term = result;
1104          for(unsigned i = start - 1; i > k; --i)
1105          {
1106             term *= ((i + 1) * y) / ((n - i) * x);
1107             result += term;
1108          }
1109          term = start_term;
1110          for(unsigned i = start + 1; i <= n; ++i)
1111          {
1112             term *= (n - i + 1) * x / (i * y);
1113             result += term;
1114          }
1115       }
1116    }
1117 
1118    return result;
1119 }
1120 
1121 
1122 //
1123 // The incomplete beta function implementation:
1124 // This is just a big bunch of spaghetti code to divide up the
1125 // input range and select the right implementation method for
1126 // each domain:
1127 //
1128 
1129 template <class T, class Policy>
1130 BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
1131 {
1132    constexpr auto function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
1133    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1134    BOOST_MATH_STD_USING // for ADL of std math functions.
1135 
1136    BOOST_MATH_INSTRUMENT_VARIABLE(a);
1137    BOOST_MATH_INSTRUMENT_VARIABLE(b);
1138    BOOST_MATH_INSTRUMENT_VARIABLE(x);
1139    BOOST_MATH_INSTRUMENT_VARIABLE(inv);
1140    BOOST_MATH_INSTRUMENT_VARIABLE(normalised);
1141 
1142    bool invert = inv;
1143    T fract;
1144    T y = 1 - x;
1145 
1146    BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
1147 
1148    if(!(boost::math::isfinite)(a))
1149       return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
1150    if(!(boost::math::isfinite)(b))
1151       return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
1152    if (!(0 <= x && x <= 1))
1153       return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
1154 
1155    if(p_derivative)
1156       *p_derivative = -1; // value not set.
1157 
1158    if(normalised)
1159    {
1160       if(a < 0)
1161          return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
1162       if(b < 0)
1163          return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
1164       // extend to a few very special cases:
1165       if(a == 0)
1166       {
1167          if(b == 0)
1168             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);
1169          if(b > 0)
1170             return static_cast<T>(inv ? 0 : 1);
1171       }
1172       else if(b == 0)
1173       {
1174          if(a > 0)
1175             return static_cast<T>(inv ? 1 : 0);
1176       }
1177    }
1178    else
1179    {
1180       if(a <= 0)
1181          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);
1182       if(b <= 0)
1183          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);
1184    }
1185 
1186    if(x == 0)
1187    {
1188       if(p_derivative)
1189       {
1190          *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
1191       }
1192       return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
1193    }
1194    if(x == 1)
1195    {
1196       if(p_derivative)
1197       {
1198          *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
1199       }
1200       return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
1201    }
1202    if((a == 0.5f) && (b == 0.5f))
1203    {
1204       // We have an arcsine distribution:
1205       if(p_derivative)
1206       {
1207          *p_derivative = 1 / (constants::pi<T>() * sqrt(y * x));
1208       }
1209       T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
1210       if(!normalised)
1211          p *= constants::pi<T>();
1212       return p;
1213    }
1214    if(a == 1)
1215    {
1216       BOOST_MATH_GPU_SAFE_SWAP(a, b);
1217       BOOST_MATH_GPU_SAFE_SWAP(x, y);
1218       invert = !invert;
1219    }
1220    if(b == 1)
1221    {
1222       //
1223       // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
1224       //
1225       if(a == 1)
1226       {
1227          if(p_derivative)
1228             *p_derivative = 1;
1229          return invert ? y : x;
1230       }
1231 
1232       if(p_derivative)
1233       {
1234          *p_derivative = a * pow(x, a - 1);
1235       }
1236       T p;  // LCOV_EXCL_LINE
1237       if(y < 0.5)
1238          p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol)));
1239       else
1240          p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a));
1241       if(!normalised)
1242          p /= a;
1243       return p;
1244    }
1245 
1246    if(BOOST_MATH_GPU_SAFE_MIN(a, b) <= 1)
1247    {
1248       if(x > 0.5)
1249       {
1250          BOOST_MATH_GPU_SAFE_SWAP(a, b);
1251          BOOST_MATH_GPU_SAFE_SWAP(x, y);
1252          invert = !invert;
1253          BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1254       }
1255       if(BOOST_MATH_GPU_SAFE_MAX(a, b) <= 1)
1256       {
1257          // Both a,b < 1:
1258          if((a >= BOOST_MATH_GPU_SAFE_MIN(T(0.2), b)) || (pow(x, a) <= 0.9))
1259          {
1260             if(!invert)
1261             {
1262                fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1263                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1264             }
1265             else
1266             {
1267                fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1268                invert = false;
1269                fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1270                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1271             }
1272          }
1273          else
1274          {
1275             BOOST_MATH_GPU_SAFE_SWAP(a, b);
1276             BOOST_MATH_GPU_SAFE_SWAP(x, y);
1277             invert = !invert;
1278             if(y >= 0.3)
1279             {
1280                if(!invert)
1281                {
1282                   fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1283                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1284                }
1285                else
1286                {
1287                   fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1288                   invert = false;
1289                   fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1290                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1291                }
1292             }
1293             else
1294             {
1295                // Sidestep on a, and then use the series representation:
1296                T prefix;  // LCOV_EXCL_LINE
1297                if(!normalised)
1298                {
1299                   prefix = rising_factorial_ratio(T(a+b), a, 20);
1300                }
1301                else
1302                {
1303                   prefix = 1;
1304                }
1305                fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1306                if(!invert)
1307                {
1308                   fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1309                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1310                }
1311                else
1312                {
1313                   fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1314                   invert = false;
1315                   fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1316                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1317                }
1318             }
1319          }
1320       }
1321       else
1322       {
1323          // One of a, b < 1 only:
1324          if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7)))
1325          {
1326             if(!invert)
1327             {
1328                fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1329                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1330             }
1331             else
1332             {
1333                fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1334                invert = false;
1335                fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1336                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1337             }
1338          }
1339          else
1340          {
1341             BOOST_MATH_GPU_SAFE_SWAP(a, b);
1342             BOOST_MATH_GPU_SAFE_SWAP(x, y);
1343             invert = !invert;
1344 
1345             if(y >= 0.3)
1346             {
1347                if(!invert)
1348                {
1349                   fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1350                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1351                }
1352                else
1353                {
1354                   fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1355                   invert = false;
1356                   fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1357                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1358                }
1359             }
1360             else if(a >= 15)
1361             {
1362                if(!invert)
1363                {
1364                   fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised);
1365                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1366                }
1367                else
1368                {
1369                   fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1370                   invert = false;
1371                   fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised);
1372                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1373                }
1374             }
1375             else
1376             {
1377                // Sidestep to improve errors:
1378                T prefix;  // LCOV_EXCL_LINE
1379                if(!normalised)
1380                {
1381                   prefix = rising_factorial_ratio(T(a+b), a, 20);
1382                }
1383                else
1384                {
1385                   prefix = 1;
1386                }
1387                fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
1388                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1389                if(!invert)
1390                {
1391                   fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1392                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1393                }
1394                else
1395                {
1396                   fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
1397                   invert = false;
1398                   fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
1399                   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1400                }
1401             }
1402          }
1403       }
1404    }
1405    else
1406    {
1407       // Both a,b >= 1:
1408       T lambda;  // LCOV_EXCL_LINE
1409       if(a < b)
1410       {
1411          lambda = a - (a + b) * x;
1412       }
1413       else
1414       {
1415          lambda = (a + b) * y - b;
1416       }
1417       if(lambda < 0)
1418       {
1419          BOOST_MATH_GPU_SAFE_SWAP(a, b);
1420          BOOST_MATH_GPU_SAFE_SWAP(x, y);
1421          invert = !invert;
1422          BOOST_MATH_INSTRUMENT_VARIABLE(invert);
1423       }
1424 
1425       if(b < 40)
1426       {
1427          if((floor(a) == a) && (floor(b) == b) && (a < static_cast<T>((boost::math::numeric_limits<int>::max)() - 100)) && (y != 1))
1428          {
1429             // relate to the binomial distribution and use a finite sum:
1430             T k = a - 1;
1431             T n = b + k;
1432             fract = binomial_ccdf(n, k, x, y, pol);
1433             if(!normalised)
1434                fract *= boost::math::beta(a, b, pol);
1435             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1436          }
1437          else if(b * x <= 0.7)
1438          {
1439             if(!invert)
1440             {
1441                fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
1442                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1443             }
1444             else
1445             {
1446                fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
1447                invert = false;
1448                fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
1449                BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1450             }
1451          }
1452          else if(a > 15)
1453          {
1454             // sidestep so we can use the series representation:
1455             int n = itrunc(T(floor(b)), pol);
1456             if(n == b)
1457                --n;
1458             T bbar = b - n;
1459             T prefix; // LCOV_EXCL_LINE
1460             if(!normalised)
1461             {
1462                prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
1463             }
1464             else
1465             {
1466                prefix = 1;
1467             }
1468             fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(nullptr));
1469             fract = beta_small_b_large_a_series(a,  bbar, x, y, fract, T(1), pol, normalised);
1470             fract /= prefix;
1471             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1472          }
1473          else if(normalised)
1474          {
1475             // The formula here for the non-normalised case is tricky to figure
1476             // out (for me!!), and requires two pochhammer calculations rather
1477             // than one, so leave it for now and only use this in the normalized case....
1478             int n = itrunc(T(floor(b)), pol);
1479             T bbar = b - n;
1480             if(bbar <= 0)
1481             {
1482                --n;
1483                bbar += 1;
1484             }
1485             fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(nullptr));
1486             fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(nullptr));
1487             if(invert)
1488                fract -= 1;  // Note this line would need changing if we ever enable this branch in non-normalized case
1489             fract = beta_small_b_large_a_series(T(a+20),  bbar, x, y, fract, T(1), pol, normalised);
1490             if(invert)
1491             {
1492                fract = -fract;
1493                invert = false;
1494             }
1495             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1496          }
1497          else
1498          {
1499             fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1500             BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1501          }
1502       }
1503       else
1504       {
1505          fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
1506          BOOST_MATH_INSTRUMENT_VARIABLE(fract);
1507       }
1508    }
1509    if(p_derivative)
1510    {
1511       if(*p_derivative < 0)
1512       {
1513          *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol);
1514       }
1515       T div = y * x;
1516 
1517       if(*p_derivative != 0)
1518       {
1519          if((tools::max_value<T>() * div < *p_derivative))
1520          {
1521             // overflow, return an arbitrarily large value:
1522             *p_derivative = tools::max_value<T>() / 2;   // LCOV_EXCL_LINE  Probably can only get here with denormalized x.
1523          }
1524          else
1525          {
1526             *p_derivative /= div;
1527          }
1528       }
1529    }
1530    return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
1531 } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)
1532 
1533 template <class T, class Policy>
1534 BOOST_MATH_GPU_ENABLED inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
1535 {
1536    return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(nullptr));
1537 }
1538 
1539 template <class T, class Policy>
1540 BOOST_MATH_GPU_ENABLED T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
1541 {
1542    constexpr auto function = "ibeta_derivative<%1%>(%1%,%1%,%1%)";
1543    //
1544    // start with the usual error checks:
1545    //
1546    if (!(boost::math::isfinite)(a))
1547       return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
1548    if (!(boost::math::isfinite)(b))
1549       return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
1550    if (!(0 <= x && x <= 1))
1551       return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
1552 
1553    if(a <= 0)
1554       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);
1555    if(b <= 0)
1556       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);
1557    //
1558    // Now the corner cases:
1559    //
1560    if(x == 0)
1561    {
1562       return (a > 1) ? 0 :
1563          (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol);
1564    }
1565    else if(x == 1)
1566    {
1567       return (b > 1) ? 0 :
1568          (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol);
1569    }
1570    //
1571    // Now the regular cases:
1572    //
1573    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1574    T y = (1 - x) * x;
1575    T f1;
1576    if (!(boost::math::isinf)(1 / y))
1577    {
1578       f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function);
1579    }
1580    else
1581    {
1582       return (a > 1) ? 0 : (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, nullptr, pol);
1583    }
1584 
1585    return f1;
1586 }
1587 //
1588 // Some forwarding functions that disambiguate the third argument type:
1589 //
1590 template <class RT1, class RT2, class Policy>
1591 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type
1592    beta(RT1 a, RT2 b, const Policy&, const boost::math::true_type*)
1593 {
1594    BOOST_FPU_EXCEPTION_GUARD
1595    typedef typename tools::promote_args<RT1, RT2>::type result_type;
1596    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1597    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1598    typedef typename policies::normalise<
1599       Policy,
1600       policies::promote_float<false>,
1601       policies::promote_double<false>,
1602       policies::discrete_quantile<>,
1603       policies::assert_undefined<> >::type forwarding_policy;
1604 
1605    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%)");
1606 }
1607 template <class RT1, class RT2, class RT3>
1608 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1609    beta(RT1 a, RT2 b, RT3 x, const boost::math::false_type*)
1610 {
1611    return boost::math::beta(a, b, x, policies::policy<>());
1612 }
1613 } // namespace detail
1614 
1615 //
1616 // The actual function entry-points now follow, these just figure out
1617 // which Lanczos approximation to use
1618 // and forward to the implementation functions:
1619 //
1620 template <class RT1, class RT2, class A>
1621 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, A>::type
1622    beta(RT1 a, RT2 b, A arg)
1623 {
1624    using tag = typename policies::is_policy<A>::type;
1625    using ReturnType = tools::promote_args_t<RT1, RT2, A>;
1626    return static_cast<ReturnType>(boost::math::detail::beta(a, b, arg, static_cast<tag*>(nullptr)));
1627 }
1628 
1629 template <class RT1, class RT2>
1630 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type
1631    beta(RT1 a, RT2 b)
1632 {
1633    return boost::math::beta(a, b, policies::policy<>());
1634 }
1635 
1636 template <class RT1, class RT2, class RT3, class Policy>
1637 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1638    beta(RT1 a, RT2 b, RT3 x, const Policy&)
1639 {
1640    BOOST_FPU_EXCEPTION_GUARD
1641    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1642    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1643    typedef typename policies::normalise<
1644       Policy,
1645       policies::promote_float<false>,
1646       policies::promote_double<false>,
1647       policies::discrete_quantile<>,
1648       policies::assert_undefined<> >::type forwarding_policy;
1649 
1650    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%)");
1651 }
1652 
1653 template <class RT1, class RT2, class RT3, class Policy>
1654 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1655    betac(RT1 a, RT2 b, RT3 x, const Policy&)
1656 {
1657    BOOST_FPU_EXCEPTION_GUARD
1658    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1659    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1660    typedef typename policies::normalise<
1661       Policy,
1662       policies::promote_float<false>,
1663       policies::promote_double<false>,
1664       policies::discrete_quantile<>,
1665       policies::assert_undefined<> >::type forwarding_policy;
1666 
1667    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%)");
1668 }
1669 template <class RT1, class RT2, class RT3>
1670 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1671    betac(RT1 a, RT2 b, RT3 x)
1672 {
1673    return boost::math::betac(a, b, x, policies::policy<>());
1674 }
1675 
1676 template <class RT1, class RT2, class RT3, class Policy>
1677 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1678    ibeta(RT1 a, RT2 b, RT3 x, const Policy&)
1679 {
1680    BOOST_FPU_EXCEPTION_GUARD
1681    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1682    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1683    typedef typename policies::normalise<
1684       Policy,
1685       policies::promote_float<false>,
1686       policies::promote_double<false>,
1687       policies::discrete_quantile<>,
1688       policies::assert_undefined<> >::type forwarding_policy;
1689 
1690    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%)");
1691 }
1692 template <class RT1, class RT2, class RT3>
1693 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1694    ibeta(RT1 a, RT2 b, RT3 x)
1695 {
1696    return boost::math::ibeta(a, b, x, policies::policy<>());
1697 }
1698 
1699 template <class RT1, class RT2, class RT3, class Policy>
1700 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1701    ibetac(RT1 a, RT2 b, RT3 x, const Policy&)
1702 {
1703    BOOST_FPU_EXCEPTION_GUARD
1704    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1705    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1706    typedef typename policies::normalise<
1707       Policy,
1708       policies::promote_float<false>,
1709       policies::promote_double<false>,
1710       policies::discrete_quantile<>,
1711       policies::assert_undefined<> >::type forwarding_policy;
1712 
1713    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%)");
1714 }
1715 template <class RT1, class RT2, class RT3>
1716 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1717    ibetac(RT1 a, RT2 b, RT3 x)
1718 {
1719    return boost::math::ibetac(a, b, x, policies::policy<>());
1720 }
1721 
1722 template <class RT1, class RT2, class RT3, class Policy>
1723 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1724    ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&)
1725 {
1726    BOOST_FPU_EXCEPTION_GUARD
1727    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1728    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1729    typedef typename policies::normalise<
1730       Policy,
1731       policies::promote_float<false>,
1732       policies::promote_double<false>,
1733       policies::discrete_quantile<>,
1734       policies::assert_undefined<> >::type forwarding_policy;
1735 
1736    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%)");
1737 }
1738 template <class RT1, class RT2, class RT3>
1739 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
1740    ibeta_derivative(RT1 a, RT2 b, RT3 x)
1741 {
1742    return boost::math::ibeta_derivative(a, b, x, policies::policy<>());
1743 }
1744 
1745 } // namespace math
1746 } // namespace boost
1747 
1748 #include <boost/math/special_functions/detail/ibeta_inverse.hpp>
1749 #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp>
1750 
1751 #endif // BOOST_MATH_SPECIAL_BETA_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!