special_functions/detail/ibeta_inverse.hpp

0001 //  Copyright John Maddock 2006.
0002 //  Copyright Paul A. Bristow 2007
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_FUNCTIONS_IBETA_INVERSE_HPP
0008 #define BOOST_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <boost/math/special_functions/beta.hpp>
0015 #include <boost/math/special_functions/erf.hpp>
0016 #include <boost/math/tools/roots.hpp>
0017 #include <boost/math/special_functions/detail/t_distribution_inv.hpp>
0018 #include <boost/math/special_functions/fpclassify.hpp>
0019 #include <boost/math/tools/precision.hpp>
0020
0021 namespace boost{ namespace math{ namespace detail{
0022
0023 //
0024 // Helper object used by root finding
0025 // code to convert eta to x.
0026 //
0027 template <class T>
0028 struct temme_root_finder
0029 {
0030    temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {
0031       BOOST_MATH_ASSERT(
0032          math::tools::epsilon<T>() <= a && !(boost::math::isinf)(a));
0033    }
0034
0035    boost::math::tuple<T, T> operator()(T x)
0036    {
0037       BOOST_MATH_STD_USING // ADL of std names
0038
0039       T y = 1 - x;
0040       T f = log(x) + a * log(y) + t;
0041       T f1 = (1 / x) - (a / (y));
0042       return boost::math::make_tuple(f, f1);
0043    }
0044 private:
0045    T t, a;
0046 };
0047 //
0048 // See:
0049 // "Asymptotic Inversion of the Incomplete Beta Function"
0050 // N.M. Temme
0051 // Journal of Computation and Applied Mathematics 41 (1992) 145-157.
0052 // Section 2.
0053 //
0054 template <class T, class Policy>
0055 T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)
0056 {
0057    BOOST_MATH_STD_USING // ADL of std names
0058
0059    const T r2 = sqrt(T(2));
0060    //
0061    // get the first approximation for eta from the inverse
0062    // error function (Eq: 2.9 and 2.10).
0063    //
0064    T eta0 = boost::math::erfc_inv(2 * z, pol);
0065    eta0 /= -sqrt(a / 2);
0066
0067    T terms[4] = { eta0 };
0068    T workspace[7];
0069    //
0070    // calculate powers:
0071    //
0072    T B = b - a;
0073    T B_2 = B * B;
0074    T B_3 = B_2 * B;
0075    //
0076    // Calculate correction terms:
0077    //
0078
0079    // See eq following 2.15:
0080    workspace[0] = -B * r2 / 2;
0081    workspace[1] = (1 - 2 * B) / 8;
0082    workspace[2] = -(B * r2 / 48);
0083    workspace[3] = T(-1) / 192;
0084    workspace[4] = -B * r2 / 3840;
0085    terms[1] = tools::evaluate_polynomial(workspace, eta0, 5);
0086    // Eq Following 2.17:
0087    workspace[0] = B * r2 * (3 * B - 2) / 12;
0088    workspace[1] = (20 * B_2 - 12 * B + 1) / 128;
0089    workspace[2] = B * r2 * (20 * B - 1) / 960;
0090    workspace[3] = (16 * B_2 + 30 * B - 15) / 4608;
0091    workspace[4] = B * r2 * (21 * B + 32) / 53760;
0092    workspace[5] = (-32 * B_2 + 63) / 368640;
0093    workspace[6] = -B * r2 * (120 * B + 17) / 25804480;
0094    terms[2] = tools::evaluate_polynomial(workspace, eta0, 7);
0095    // Eq Following 2.17:
0096    workspace[0] = B * r2 * (-75 * B_2 + 80 * B - 16) / 480;
0097    workspace[1] = (-1080 * B_3 + 868 * B_2 - 90 * B - 45) / 9216;
0098    workspace[2] = B * r2 * (-1190 * B_2 + 84 * B + 373) / 53760;
0099    workspace[3] = (-2240 * B_3 - 2508 * B_2 + 2100 * B - 165) / 368640;
0100    terms[3] = tools::evaluate_polynomial(workspace, eta0, 4);
0101    //
0102    // Bring them together to get a final estimate for eta:
0103    //
0104    T eta = tools::evaluate_polynomial(terms, T(1/a), 4);
0105    //
0106    // now we need to convert eta to x, by solving the appropriate
0107    // quadratic equation:
0108    //
0109    T eta_2 = eta * eta;
0110    T c = -exp(-eta_2 / 2);
0111    T x;
0112    if(eta_2 == 0)
0113       x = static_cast<T>(0.5f);
0114    else
0115       x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;
0116    //
0117    // These are post-conditions of the method, but the addition above
0118    // may result in us being out by 1ulp either side of the boundary,
0119    // so just check that we're in bounds and adjust as needed.
0120    // See https://github.com/boostorg/math/issues/961
0121    //
0122    if (x < 0)
0123       x = 0;
0124    else if (x > 1)
0125       x = 1;
0126
0127    BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0);
0128 #ifdef BOOST_INSTRUMENT
0129    std::cout << "Estimating x with Temme method 1: " << x << std::endl;
0130 #endif
0131    return x;
0132 }
0133 //
0134 // See:
0135 // "Asymptotic Inversion of the Incomplete Beta Function"
0136 // N.M. Temme
0137 // Journal of Computation and Applied Mathematics 41 (1992) 145-157.
0138 // Section 3.
0139 //
0140 template <class T, class Policy>
0141 T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r, T theta, const Policy& pol)
0142 {
0143    BOOST_MATH_STD_USING // ADL of std names
0144
0145    //
0146    // Get first estimate for eta, see Eq 3.9 and 3.10,
0147    // but note there is a typo in Eq 3.10:
0148    //
0149    T eta0 = boost::math::erfc_inv(2 * z, pol);
0150    eta0 /= -sqrt(r / 2);
0151
0152    T s = sin(theta);
0153    T c = cos(theta);
0154    //
0155    // Now we need to perturb eta0 to get eta, which we do by
0156    // evaluating the polynomial in 1/r at the bottom of page 151,
0157    // to do this we first need the error terms e1, e2 e3
0158    // which we'll fill into the array "terms".  Since these
0159    // terms are themselves polynomials, we'll need another
0160    // array "workspace" to calculate those...
0161    //
0162    T terms[4] = { eta0 };
0163    T workspace[6];
0164    //
0165    // some powers of sin(theta)cos(theta) that we'll need later:
0166    //
0167    T sc = s * c;
0168    T sc_2 = sc * sc;
0169    T sc_3 = sc_2 * sc;
0170    T sc_4 = sc_2 * sc_2;
0171    T sc_5 = sc_2 * sc_3;
0172    T sc_6 = sc_3 * sc_3;
0173    T sc_7 = sc_4 * sc_3;
0174    //
0175    // Calculate e1 and put it in terms[1], see the middle of page 151:
0176    //
0177    workspace[0] = (2 * s * s - 1) / (3 * s * c);
0178    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co1[] = { -1, -5, 5 };
0179    workspace[1] = -tools::evaluate_even_polynomial(co1, s, 3) / (36 * sc_2);
0180    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co2[] = { 1, 21, -69, 46 };
0181    workspace[2] = tools::evaluate_even_polynomial(co2, s, 4) / (1620 * sc_3);
0182    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co3[] = { 7, -2, 33, -62, 31 };
0183    workspace[3] = -tools::evaluate_even_polynomial(co3, s, 5) / (6480 * sc_4);
0184    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co4[] = { 25, -52, -17, 88, -115, 46 };
0185    workspace[4] = tools::evaluate_even_polynomial(co4, s, 6) / (90720 * sc_5);
0186    terms[1] = tools::evaluate_polynomial(workspace, eta0, 5);
0187    //
0188    // Now evaluate e2 and put it in terms[2]:
0189    //
0190    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co5[] = { 7, 12, -78, 52 };
0191    workspace[0] = -tools::evaluate_even_polynomial(co5, s, 4) / (405 * sc_3);
0192    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co6[] = { -7, 2, 183, -370, 185 };
0193    workspace[1] = tools::evaluate_even_polynomial(co6, s, 5) / (2592 * sc_4);
0194    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co7[] = { -533, 776, -1835, 10240, -13525, 5410 };
0195    workspace[2] = -tools::evaluate_even_polynomial(co7, s, 6) / (204120 * sc_5);
0196    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co8[] = { -1579, 3747, -3372, -15821, 45588, -45213, 15071 };
0197    workspace[3] = -tools::evaluate_even_polynomial(co8, s, 7) / (2099520 * sc_6);
0198    terms[2] = tools::evaluate_polynomial(workspace, eta0, 4);
0199    //
0200    // And e3, and put it in terms[3]:
0201    //
0202    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co9[] = {449, -1259, -769, 6686, -9260, 3704 };
0203    workspace[0] = tools::evaluate_even_polynomial(co9, s, 6) / (102060 * sc_5);
0204    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co10[] = { 63149, -151557, 140052, -727469, 2239932, -2251437, 750479 };
0205    workspace[1] = -tools::evaluate_even_polynomial(co10, s, 7) / (20995200 * sc_6);
0206    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co11[] = { 29233, -78755, 105222, 146879, -1602610, 3195183, -2554139, 729754 };
0207    workspace[2] = tools::evaluate_even_polynomial(co11, s, 8) / (36741600 * sc_7);
0208    terms[3] = tools::evaluate_polynomial(workspace, eta0, 3);
0209    //
0210    // Bring the correction terms together to evaluate eta,
0211    // this is the last equation on page 151:
0212    //
0213    T eta = tools::evaluate_polynomial(terms, T(1/r), 4);
0214    //
0215    // Now that we have eta we need to back solve for x,
0216    // we seek the value of x that gives eta in Eq 3.2.
0217    // The two methods used are described in section 5.
0218    //
0219    // Begin by defining a few variables we'll need later:
0220    //
0221    T x;
0222    T s_2 = s * s;
0223    T c_2 = c * c;
0224    T alpha = c / s;
0225    alpha *= alpha;
0226    T lu = (-(eta * eta) / (2 * s_2) + log(s_2) + c_2 * log(c_2) / s_2);
0227    //
0228    // Temme doesn't specify what value to switch on here,
0229    // but this seems to work pretty well:
0230    //
0231    if(fabs(eta) < 0.7)
0232    {
0233       //
0234       // Small eta use the expansion Temme gives in the second equation
0235       // of section 5, it's a polynomial in eta:
0236       //
0237       workspace[0] = s * s;
0238       workspace[1] = s * c;
0239       workspace[2] = (1 - 2 * workspace[0]) / 3;
0240       static const BOOST_MATH_INT_TABLE_TYPE(T, int) co12[] = { 1, -13, 13 };
0241       workspace[3] = tools::evaluate_polynomial(co12, workspace[0], 3) / (36 * s * c);
0242       static const BOOST_MATH_INT_TABLE_TYPE(T, int) co13[] = { 1, 21, -69, 46 };
0243       workspace[4] = tools::evaluate_polynomial(co13, workspace[0], 4) / (270 * workspace[0] * c * c);
0244       x = tools::evaluate_polynomial(workspace, eta, 5);
0245 #ifdef BOOST_INSTRUMENT
0246       std::cout << "Estimating x with Temme method 2 (small eta): " << x << std::endl;
0247 #endif
0248    }
0249    else
0250    {
0251       //
0252       // If eta is large we need to solve Eq 3.2 more directly,
0253       // begin by getting an initial approximation for x from
0254       // the last equation on page 155, this is a polynomial in u:
0255       //
0256       T u = exp(lu);
0257       workspace[0] = u;
0258       workspace[1] = alpha;
0259       workspace[2] = 0;
0260       workspace[3] = 3 * alpha * (3 * alpha + 1) / 6;
0261       workspace[4] = 4 * alpha * (4 * alpha + 1) * (4 * alpha + 2) / 24;
0262       workspace[5] = 5 * alpha * (5 * alpha + 1) * (5 * alpha + 2) * (5 * alpha + 3) / 120;
0263       x = tools::evaluate_polynomial(workspace, u, 6);
0264       //
0265       // At this point we may or may not have the right answer, Eq-3.2 has
0266       // two solutions for x for any given eta, however the mapping in 3.2
0267       // is 1:1 with the sign of eta and x-sin^2(theta) being the same.
0268       // So we can check if we have the right root of 3.2, and if not
0269       // switch x for 1-x.  This transformation is motivated by the fact
0270       // that the distribution is *almost* symmetric so 1-x will be in the right
0271       // ball park for the solution:
0272       //
0273       if((x - s_2) * eta < 0)
0274          x = 1 - x;
0275 #ifdef BOOST_INSTRUMENT
0276       std::cout << "Estimating x with Temme method 2 (large eta): " << x << std::endl;
0277 #endif
0278    }
0279    //
0280    // The final step is a few Newton-Raphson iterations to
0281    // clean up our approximation for x, this is pretty cheap
0282    // in general, and very cheap compared to an incomplete beta
0283    // evaluation.  The limits set on x come from the observation
0284    // that the sign of eta and x-sin^2(theta) are the same.
0285    //
0286    T lower, upper;
0287    if(eta < 0)
0288    {
0289       lower = 0;
0290       upper = s_2;
0291    }
0292    else
0293    {
0294       lower = s_2;
0295       upper = 1;
0296    }
0297    //
0298    // If our initial approximation is out of bounds then bisect:
0299    //
0300    if((x < lower) || (x > upper))
0301       x = (lower+upper) / 2;
0302    //
0303    // And iterate:
0304    //
0305    x = tools::newton_raphson_iterate(
0306       temme_root_finder<T>(-lu, alpha), x, lower, upper, policies::digits<T, Policy>() / 2);
0307
0308    return x;
0309 }
0310 //
0311 // See:
0312 // "Asymptotic Inversion of the Incomplete Beta Function"
0313 // N.M. Temme
0314 // Journal of Computation and Applied Mathematics 41 (1992) 145-157.
0315 // Section 4.
0316 //
0317 template <class T, class Policy>
0318 T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
0319 {
0320    BOOST_MATH_STD_USING // ADL of std names
0321
0322    //
0323    // Begin by getting an initial approximation for the quantity
0324    // eta from the dominant part of the incomplete beta:
0325    //
0326    T eta0;
0327    if(p < q)
0328       eta0 = boost::math::gamma_q_inv(b, p, pol);
0329    else
0330       eta0 = boost::math::gamma_p_inv(b, q, pol);
0331    eta0 /= a;
0332    //
0333    // Define the variables and powers we'll need later on:
0334    //
0335    T mu = b / a;
0336    T w = sqrt(1 + mu);
0337    T w_2 = w * w;
0338    T w_3 = w_2 * w;
0339    T w_4 = w_2 * w_2;
0340    T w_5 = w_3 * w_2;
0341    T w_6 = w_3 * w_3;
0342    T w_7 = w_4 * w_3;
0343    T w_8 = w_4 * w_4;
0344    T w_9 = w_5 * w_4;
0345    T w_10 = w_5 * w_5;
0346    T d = eta0 - mu;
0347    T d_2 = d * d;
0348    T d_3 = d_2 * d;
0349    T d_4 = d_2 * d_2;
0350    T w1 = w + 1;
0351    T w1_2 = w1 * w1;
0352    T w1_3 = w1 * w1_2;
0353    T w1_4 = w1_2 * w1_2;
0354    //
0355    // Now we need to compute the perturbation error terms that
0356    // convert eta0 to eta, these are all polynomials of polynomials.
0357    // Probably these should be re-written to use tabulated data
0358    // (see examples above), but it's less of a win in this case as we
0359    // need to calculate the individual powers for the denominator terms
0360    // anyway, so we might as well use them for the numerator-polynomials
0361    // as well....
0362    //
0363    // Refer to p154-p155 for the details of these expansions:
0364    //
0365    T e1 = (w + 2) * (w - 1) / (3 * w);
0366    e1 += (w_3 + 9 * w_2 + 21 * w + 5) * d / (36 * w_2 * w1);
0367    e1 -= (w_4 - 13 * w_3 + 69 * w_2 + 167 * w + 46) * d_2 / (1620 * w1_2 * w_3);
0368    e1 -= (7 * w_5 + 21 * w_4 + 70 * w_3 + 26 * w_2 - 93 * w - 31) * d_3 / (6480 * w1_3 * w_4);
0369    e1 -= (75 * w_6 + 202 * w_5 + 188 * w_4 - 888 * w_3 - 1345 * w_2 + 118 * w + 138) * d_4 / (272160 * w1_4 * w_5);
0370
0371    T e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1) / (1620 * w1 * w_3);
0372    e2 -= (35 * w_6 - 154 * w_5 - 623 * w_4 - 1636 * w_3 - 3983 * w_2 - 3514 * w - 925) * d / (12960 * w1_2 * w_4);
0373    e2 -= (2132 * w_7 + 7915 * w_6 + 16821 * w_5 + 35066 * w_4 + 87490 * w_3 + 141183 * w_2 + 95993 * w + 21640) * d_2  / (816480 * w_5 * w1_3);
0374    e2 -= (11053 * w_8 + 53308 * w_7 + 117010 * w_6 + 163924 * w_5 + 116188 * w_4 - 258428 * w_3 - 677042 * w_2 - 481940 * w - 105497) * d_3 / (T(14696640) * w1_4 * w_6);
0375
0376    T e3 = -((3592 * w_7 + 8375 * w_6 - 1323 * w_5 - 29198 * w_4 - 89578 * w_3 - 154413 * w_2 - 116063 * w - 29632) * (w - 1)) / (816480 * w_5 * w1_2);
0377    e3 -= (442043 * w_9 + T(2054169) * w_8 + T(3803094) * w_7 + T(3470754) * w_6 + T(2141568) * w_5 - T(2393568) * w_4 - T(19904934) * w_3 - T(34714674) * w_2 - T(23128299) * w - T(5253353)) * d / (T(146966400) * w_6 * w1_3);
0378    e3 -= (116932 * w_10 + 819281 * w_9 + T(2378172) * w_8 + T(4341330) * w_7 + T(6806004) * w_6 + T(10622748) * w_5 + T(18739500) * w_4 + T(30651894) * w_3 + T(30869976) * w_2 + T(15431867) * w + T(2919016)) * d_2 / (T(146966400) * w1_4 * w_7);
0379    //
0380    // Combine eta0 and the error terms to compute eta (Second equation p155):
0381    //
0382    T eta = eta0 + e1 / a + e2 / (a * a) + e3 / (a * a * a);
0383    //
0384    // Now we need to solve Eq 4.2 to obtain x.  For any given value of
0385    // eta there are two solutions to this equation, and since the distribution
0386    // may be very skewed, these are not related by x ~ 1-x we used when
0387    // implementing section 3 above.  However we know that:
0388    //
0389    //  cross < x <= 1       ; iff eta < mu
0390    //          x == cross   ; iff eta == mu
0391    //     0 <= x < cross    ; iff eta > mu
0392    //
0393    // Where cross == 1 / (1 + mu)
0394    // Many thanks to Prof Temme for clarifying this point.
0395    //
0396    // Therefore we'll just jump straight into Newton iterations
0397    // to solve Eq 4.2 using these bounds, and simple bisection
0398    // as the first guess, in practice this converges pretty quickly
0399    // and we only need a few digits correct anyway:
0400    //
0401    if(eta <= 0)
0402       eta = tools::min_value<T>();
0403    T u = eta - mu * log(eta) + (1 + mu) * log(1 + mu) - mu;
0404    T cross = 1 / (1 + mu);
0405    T lower = eta < mu ? cross : 0;
0406    T upper = eta < mu ? 1 : cross;
0407    T x = (lower + upper) / 2;
0408
0409    // Early exit for cases with numerical precision issues.
0410    if (cross == 0 || cross == 1) { return cross; }
0411
0412    x = tools::newton_raphson_iterate(
0413       temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2);
0414 #ifdef BOOST_INSTRUMENT
0415    std::cout << "Estimating x with Temme method 3: " << x << std::endl;
0416 #endif
0417    return x;
0418 }
0419
0420 template <class T, class Policy>
0421 struct ibeta_roots
0422 {
0423    ibeta_roots(T _a, T _b, T t, bool inv = false)
0424       : a(_a), b(_b), target(t), invert(inv) {}
0425
0426    boost::math::tuple<T, T, T> operator()(T x)
0427    {
0428       BOOST_MATH_STD_USING // ADL of std names
0429
0430       BOOST_FPU_EXCEPTION_GUARD
0431
0432       T f1;
0433       T y = 1 - x;
0434       T f = ibeta_imp(a, b, x, Policy(), invert, true, &f1) - target;
0435       if(invert)
0436          f1 = -f1;
0437       if(y == 0)
0438          y = tools::min_value<T>() * 64;
0439       if(x == 0)
0440          x = tools::min_value<T>() * 64;
0441
0442       T f2 = f1 * (-y * a + (b - 2) * x + 1);
0443       if(fabs(f2) < y * x * tools::max_value<T>())
0444          f2 /= (y * x);
0445       if(invert)
0446          f2 = -f2;
0447
0448       // make sure we don't have a zero derivative:
0449       if(f1 == 0)
0450          f1 = (invert ? -1 : 1) * tools::min_value<T>() * 64;
0451
0452       return boost::math::make_tuple(f, f1, f2);
0453    }
0454 private:
0455    T a, b, target;
0456    bool invert;
0457 };
0458
0459 template <class T, class Policy>
0460 T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
0461 {
0462    BOOST_MATH_STD_USING  // For ADL of math functions.
0463
0464    //
0465    // The flag invert is set to true if we swap a for b and p for q,
0466    // in which case the result has to be subtracted from 1:
0467    //
0468    bool invert = false;
0469    //
0470    // Handle trivial cases first:
0471    //
0472    if(q == 0)
0473    {
0474       if(py) *py = 0;
0475       return 1;
0476    }
0477    else if(p == 0)
0478    {
0479       if(py) *py = 1;
0480       return 0;
0481    }
0482    else if(a == 1)
0483    {
0484       if(b == 1)
0485       {
0486          if(py) *py = 1 - p;
0487          return p;
0488       }
0489       // Change things around so we can handle as b == 1 special case below:
0490       std::swap(a, b);
0491       std::swap(p, q);
0492       invert = true;
0493    }
0494    //
0495    // Depending upon which approximation method we use, we may end up
0496    // calculating either x or y initially (where y = 1-x):
0497    //
0498    T x = 0; // Set to a safe zero to avoid a
0499    // MSVC 2005 warning C4701: potentially uninitialized local variable 'x' used
0500    // But code inspection appears to ensure that x IS assigned whatever the code path.
0501    T y;
0502
0503    // For some of the methods we can put tighter bounds
0504    // on the result than simply [0,1]:
0505    //
0506    T lower = 0;
0507    T upper = 1;
0508    //
0509    // Student's T with b = 0.5 gets handled as a special case, swap
0510    // around if the arguments are in the "wrong" order:
0511    //
0512    if(a == 0.5f)
0513    {
0514       if(b == 0.5f)
0515       {
0516          x = sin(p * constants::half_pi<T>());
0517          x *= x;
0518          if(py)
0519          {
0520             *py = sin(q * constants::half_pi<T>());
0521             *py *= *py;
0522          }
0523          return x;
0524       }
0525       else if(b > 0.5f)
0526       {
0527          std::swap(a, b);
0528          std::swap(p, q);
0529          invert = !invert;
0530       }
0531    }
0532    //
0533    // Select calculation method for the initial estimate:
0534    //
0535    if((b == 0.5f) && (a >= 0.5f) && (p != 1))
0536    {
0537       //
0538       // We have a Student's T distribution:
0539       x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol);
0540    }
0541    else if(b == 1)
0542    {
0543       if(p < q)
0544       {
0545          if(a > 1)
0546          {
0547             x = pow(p, 1 / a);
0548             y = -boost::math::expm1(log(p) / a, pol);
0549          }
0550          else
0551          {
0552             x = pow(p, 1 / a);
0553             y = 1 - x;
0554          }
0555       }
0556       else
0557       {
0558          x = exp(boost::math::log1p(-q, pol) / a);
0559          y = -boost::math::expm1(boost::math::log1p(-q, pol) / a, pol);
0560       }
0561       if(invert)
0562          std::swap(x, y);
0563       if(py)
0564          *py = y;
0565       return x;
0566    }
0567    else if(a + b > 5)
0568    {
0569       //
0570       // When a+b is large then we can use one of Prof Temme's
0571       // asymptotic expansions, begin by swapping things around
0572       // so that p < 0.5, we do this to avoid cancellations errors
0573       // when p is large.
0574       //
0575       if(p > 0.5)
0576       {
0577          std::swap(a, b);
0578          std::swap(p, q);
0579          invert = !invert;
0580       }
0581       T minv = (std::min)(a, b);
0582       T maxv = (std::max)(a, b);
0583       if((sqrt(minv) > (maxv - minv)) && (minv > 5))
0584       {
0585          //
0586          // When a and b differ by a small amount
0587          // the curve is quite symmetrical and we can use an error
0588          // function to approximate the inverse. This is the cheapest
0589          // of the three Temme expansions, and the calculated value
0590          // for x will never be much larger than p, so we don't have
0591          // to worry about cancellation as long as p is small.
0592          //
0593          x = temme_method_1_ibeta_inverse(a, b, p, pol);
0594          y = 1 - x;
0595       }
0596       else
0597       {
0598          T r = a + b;
0599          T theta = asin(sqrt(a / r));
0600          T lambda = minv / r;
0601          if((lambda >= 0.2) && (lambda <= 0.8) && (r >= 10))
0602          {
0603             //
0604             // The second error function case is the next cheapest
0605             // to use, it brakes down when the result is likely to be
0606             // very small, if a+b is also small, but we can use a
0607             // cheaper expansion there in any case.  As before x won't
0608             // be much larger than p, so as long as p is small we should
0609             // be free of cancellation error.
0610             //
0611             T ppa = pow(p, 1/a);
0612             if((ppa < 0.0025) && (a + b < 200))
0613             {
0614                x = ppa * pow(a * boost::math::beta(a, b, pol), 1/a);
0615             }
0616             else
0617                x = temme_method_2_ibeta_inverse(a, b, p, r, theta, pol);
0618             y = 1 - x;
0619          }
0620          else
0621          {
0622             //
0623             // If we get here then a and b are very different in magnitude
0624             // and we need to use the third of Temme's methods which
0625             // involves inverting the incomplete gamma.  This is much more
0626             // expensive than the other methods.  We also can only use this
0627             // method when a > b, which can lead to cancellation errors
0628             // if we really want y (as we will when x is close to 1), so
0629             // a different expansion is used in that case.
0630             //
0631             if(a < b)
0632             {
0633                std::swap(a, b);
0634                std::swap(p, q);
0635                invert = !invert;
0636             }
0637             //
0638             // Try and compute the easy way first:
0639             //
0640             T bet = 0;
0641             if (b < 2)
0642             {
0643 #ifndef BOOST_MATH_NO_EXCEPTIONS
0644                try
0645 #endif
0646                {
0647                   bet = boost::math::beta(a, b, pol);
0648
0649                   typedef typename Policy::overflow_error_type overflow_type;
0650
0651                   BOOST_MATH_IF_CONSTEXPR(overflow_type::value != boost::math::policies::throw_on_error)
0652                      if(bet > tools::max_value<T>())
0653                         bet = tools::max_value<T>();
0654                }
0655 #ifndef BOOST_MATH_NO_EXCEPTIONS
0656                catch (const std::overflow_error&)
0657                {
0658                   bet = tools::max_value<T>();
0659                }
0660 #endif
0661             }
0662             if(bet != 0)
0663             {
0664                y = pow(b * q * bet, 1/b);
0665                x = 1 - y;
0666             }
0667             else
0668                y = 1;
0669             if(y > 1e-5)
0670             {
0671                x = temme_method_3_ibeta_inverse(a, b, p, q, pol);
0672                y = 1 - x;
0673             }
0674          }
0675       }
0676    }
0677    else if((a < 1) && (b < 1))
0678    {
0679       //
0680       // Both a and b less than 1,
0681       // there is a point of inflection at xs:
0682       //
0683       T xs = (1 - a) / (2 - a - b);
0684       //
0685       // Now we need to ensure that we start our iteration from the
0686       // right side of the inflection point:
0687       //
0688       T fs = boost::math::ibeta(a, b, xs, pol) - p;
0689       if(fabs(fs) / p < tools::epsilon<T>() * 3)
0690       {
0691          // The result is at the point of inflection, best just return it:
0692          *py = invert ? xs : 1 - xs;
0693          return invert ? 1-xs : xs;
0694       }
0695       if(fs < 0)
0696       {
0697          std::swap(a, b);
0698          std::swap(p, q);
0699          invert = !invert;
0700          xs = 1 - xs;
0701       }
0702       if ((a < tools::min_value<T>()) && (b > tools::min_value<T>()))
0703       {
0704          if (py)
0705          {
0706             *py = invert ? 0 : 1;
0707          }
0708          return invert ? 1 : 0; // nothing interesting going on here.
0709       }
0710       //
0711       // The call to beta may overflow, plus the alternative using lgamma may do the same
0712       // if T is a type where 1/T is infinite for small values (denorms for example).
0713       //
0714       T bet = 0;
0715       T xg;
0716       bool overflow = false;
0717 #ifndef BOOST_MATH_NO_EXCEPTIONS
0718       try {
0719 #endif
0720          bet = boost::math::beta(a, b, pol);
0721 #ifndef BOOST_MATH_NO_EXCEPTIONS
0722       }
0723       catch (const std::runtime_error&)
0724       {
0725          overflow = true;
0726       }
0727 #endif
0728       if (overflow || !(boost::math::isfinite)(bet))
0729       {
0730          xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a);
0731          if (xg > 2 / tools::epsilon<T>())
0732             xg = 2 / tools::epsilon<T>();
0733       }
0734       else
0735          xg = pow(a * p * bet, 1/a);
0736       x = xg / (1 + xg);
0737       y = 1 / (1 + xg);
0738       //
0739       // And finally we know that our result is below the inflection
0740       // point, so set an upper limit on our search:
0741       //
0742       if(x > xs)
0743          x = xs;
0744       upper = xs;
0745    }
0746    else if((a > 1) && (b > 1))
0747    {
0748       //
0749       // Small a and b, both greater than 1,
0750       // there is a point of inflection at xs,
0751       // and it's complement is xs2, we must always
0752       // start our iteration from the right side of the
0753       // point of inflection.
0754       //
0755       T xs = (a - 1) / (a + b - 2);
0756       T xs2 = (b - 1) / (a + b - 2);
0757       T ps = boost::math::ibeta(a, b, xs, pol) - p;
0758
0759       if(ps < 0)
0760       {
0761          std::swap(a, b);
0762          std::swap(p, q);
0763          std::swap(xs, xs2);
0764          invert = !invert;
0765       }
0766       //
0767       // Estimate x and y, using expm1 to get a good estimate
0768       // for y when it's very small:
0769       //
0770       T lx = log(p * a * boost::math::beta(a, b, pol)) / a;
0771       x = exp(lx);
0772       y = x < 0.9 ? T(1 - x) : (T)(-boost::math::expm1(lx, pol));
0773
0774       if((b < a) && (x < 0.2))
0775       {
0776          //
0777          // Under a limited range of circumstances we can improve
0778          // our estimate for x, frankly it's clear if this has much effect!
0779          //
0780          T ap1 = a - 1;
0781          T bm1 = b - 1;
0782          T a_2 = a * a;
0783          T a_3 = a * a_2;
0784          T b_2 = b * b;
0785          T terms[5] = { 0, 1 };
0786          terms[2] = bm1 / ap1;
0787          ap1 *= ap1;
0788          terms[3] = bm1 * (3 * a * b + 5 * b + a_2 - a - 4) / (2 * (a + 2) * ap1);
0789          ap1 *= (a + 1);
0790          terms[4] = bm1 * (33 * a * b_2 + 31 * b_2 + 8 * a_2 * b_2 - 30 * a * b - 47 * b + 11 * a_2 * b + 6 * a_3 * b + 18 + 4 * a - a_3 + a_2 * a_2 - 10 * a_2)
0791                     / (3 * (a + 3) * (a + 2) * ap1);
0792          x = tools::evaluate_polynomial(terms, x, 5);
0793       }
0794       //
0795       // And finally we know that our result is below the inflection
0796       // point, so set an upper limit on our search:
0797       //
0798       if(x > xs)
0799          x = xs;
0800       upper = xs;
0801    }
0802    else /*if((a <= 1) != (b <= 1))*/
0803    {
0804       //
0805       // If all else fails we get here, only one of a and b
0806       // is above 1, and a+b is small.  Start by swapping
0807       // things around so that we have a concave curve with b > a
0808       // and no points of inflection in [0,1].  As long as we expect
0809       // x to be small then we can use the simple (and cheap) power
0810       // term to estimate x, but when we expect x to be large then
0811       // this greatly underestimates x and leaves us trying to
0812       // iterate "round the corner" which may take almost forever...
0813       //
0814       // We could use Temme's inverse gamma function case in that case,
0815       // this works really rather well (albeit expensively) even though
0816       // strictly speaking we're outside it's defined range.
0817       //
0818       // However it's expensive to compute, and an alternative approach
0819       // which models the curve as a distorted quarter circle is much
0820       // cheaper to compute, and still keeps the number of iterations
0821       // required down to a reasonable level.  With thanks to Prof Temme
0822       // for this suggestion.
0823       //
0824       if(b < a)
0825       {
0826          std::swap(a, b);
0827          std::swap(p, q);
0828          invert = !invert;
0829       }
0830       if (a < tools::min_value<T>())
0831       {
0832          // Avoid spurious overflows for denorms:
0833          if (p < 1)
0834          {
0835             x = 1;
0836             y = 0;
0837          }
0838          else
0839          {
0840             x = 0;
0841             y = 1;
0842          }
0843       }
0844       else if(pow(p, 1/a) < 0.5)
0845       {
0846 #ifndef BOOST_MATH_NO_EXCEPTIONS
0847          try
0848          {
0849 #endif
0850             x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
0851             if ((x > 1) || !(boost::math::isfinite)(x))
0852                x = 1;
0853 #ifndef BOOST_MATH_NO_EXCEPTIONS
0854          }
0855          catch (const std::overflow_error&)
0856          {
0857             x = 1;
0858          }
0859 #endif
0860          if(x == 0)
0861             x = boost::math::tools::min_value<T>();
0862          y = 1 - x;
0863       }
0864       else /*if(pow(q, 1/b) < 0.1)*/
0865       {
0866          // model a distorted quarter circle:
0867 #ifndef BOOST_MATH_NO_EXCEPTIONS
0868          try
0869          {
0870 #endif
0871             y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
0872             if ((y > 1) || !(boost::math::isfinite)(y))
0873                y = 1;
0874 #ifndef BOOST_MATH_NO_EXCEPTIONS
0875          }
0876          catch (const std::overflow_error&)
0877          {
0878             y = 1;
0879          }
0880 #endif
0881          if(y == 0)
0882             y = boost::math::tools::min_value<T>();
0883          x = 1 - y;
0884       }
0885    }
0886
0887    //
0888    // Now we have a guess for x (and for y) we can set things up for
0889    // iteration.  If x > 0.5 it pays to swap things round:
0890    //
0891    if(x > 0.5)
0892    {
0893       std::swap(a, b);
0894       std::swap(p, q);
0895       std::swap(x, y);
0896       invert = !invert;
0897       T l = 1 - upper;
0898       T u = 1 - lower;
0899       lower = l;
0900       upper = u;
0901    }
0902    //
0903    // lower bound for our search:
0904    //
0905    // We're not interested in denormalised answers as these tend to
0906    // these tend to take up lots of iterations, given that we can't get
0907    // accurate derivatives in this area (they tend to be infinite).
0908    //
0909    if(lower == 0)
0910    {
0911       if(invert && (py == 0))
0912       {
0913          //
0914          // We're not interested in answers smaller than machine epsilon:
0915          //
0916          lower = boost::math::tools::epsilon<T>();
0917          if(x < lower)
0918             x = lower;
0919       }
0920       else
0921          lower = boost::math::tools::min_value<T>();
0922       if(x < lower)
0923          x = lower;
0924    }
0925    std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0926    std::uintmax_t max_iter_used = 0;
0927    //
0928    // Figure out how many digits to iterate towards:
0929    //
0930    int digits = boost::math::policies::digits<T, Policy>() / 2;
0931    if((x < 1e-50) && ((a < 1) || (b < 1)))
0932    {
0933       //
0934       // If we're in a region where the first derivative is very
0935       // large, then we have to take care that the root-finder
0936       // doesn't terminate prematurely.  We'll bump the precision
0937       // up to avoid this, but we have to take care not to set the
0938       // precision too high or the last few iterations will just
0939       // thrash around and convergence may be slow in this case.
0940       // Try 3/4 of machine epsilon:
0941       //
0942       digits *= 3;
0943       digits /= 2;
0944    }
0945    //
0946    // Now iterate, we can use either p or q as the target here
0947    // depending on which is smaller:
0948    //
0949    x = boost::math::tools::halley_iterate(
0950       boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
0951    policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol);
0952    //
0953    // We don't really want these asserts here, but they are useful for sanity
0954    // checking that we have the limits right, uncomment if you suspect bugs *only*.
0955    //
0956    //BOOST_MATH_ASSERT(x != upper);
0957    //BOOST_MATH_ASSERT((x != lower) || (x == boost::math::tools::min_value<T>()) || (x == boost::math::tools::epsilon<T>()));
0958    //
0959    // Tidy up, if we "lower" was too high then zero is the best answer we have:
0960    //
0961    if(x == lower)
0962       x = 0;
0963    if(py)
0964       *py = invert ? x : 1 - x;
0965    return invert ? 1-x : x;
0966 }
0967
0968 } // namespace detail
0969
0970 template <class T1, class T2, class T3, class T4, class Policy>
0971 inline typename tools::promote_args<T1, T2, T3, T4>::type
0972    ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy& pol)
0973 {
0974    static const char* function = "boost::math::ibeta_inv<%1%>(%1%,%1%,%1%)";
0975    BOOST_FPU_EXCEPTION_GUARD
0976    typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type;
0977    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0978    typedef typename policies::normalise<
0979       Policy,
0980       policies::promote_float<false>,
0981       policies::promote_double<false>,
0982       policies::discrete_quantile<>,
0983       policies::assert_undefined<> >::type forwarding_policy;
0984
0985    if(a <= 0)
0986       return policies::raise_domain_error<result_type>(function, "The argument a to the incomplete beta function inverse must be greater than zero (got a=%1%).", a, pol);
0987    if(b <= 0)
0988       return policies::raise_domain_error<result_type>(function, "The argument b to the incomplete beta function inverse must be greater than zero (got b=%1%).", b, pol);
0989    if((p < 0) || (p > 1))
0990       return policies::raise_domain_error<result_type>(function, "Argument p outside the range [0,1] in the incomplete beta function inverse (got p=%1%).", p, pol);
0991
0992    value_type rx, ry;
0993
0994    rx = detail::ibeta_inv_imp(
0995          static_cast<value_type>(a),
0996          static_cast<value_type>(b),
0997          static_cast<value_type>(p),
0998          static_cast<value_type>(1 - p),
0999          forwarding_policy(), &ry);
1000
1001    if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function);
1002    return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function);
1003 }
1004
1005 template <class T1, class T2, class T3, class T4>
1006 inline typename tools::promote_args<T1, T2, T3, T4>::type
1007    ibeta_inv(T1 a, T2 b, T3 p, T4* py)
1008 {
1009    return ibeta_inv(a, b, p, py, policies::policy<>());
1010 }
1011
1012 template <class T1, class T2, class T3>
1013 inline typename tools::promote_args<T1, T2, T3>::type
1014    ibeta_inv(T1 a, T2 b, T3 p)
1015 {
1016    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
1017    return ibeta_inv(a, b, p, static_cast<result_type*>(nullptr), policies::policy<>());
1018 }
1019
1020 template <class T1, class T2, class T3, class Policy>
1021 inline typename tools::promote_args<T1, T2, T3>::type
1022    ibeta_inv(T1 a, T2 b, T3 p, const Policy& pol)
1023 {
1024    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
1025    return ibeta_inv(a, b, p, static_cast<result_type*>(nullptr), pol);
1026 }
1027
1028 template <class T1, class T2, class T3, class T4, class Policy>
1029 inline typename tools::promote_args<T1, T2, T3, T4>::type
1030    ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy& pol)
1031 {
1032    static const char* function = "boost::math::ibetac_inv<%1%>(%1%,%1%,%1%)";
1033    BOOST_FPU_EXCEPTION_GUARD
1034    typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type;
1035    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1036    typedef typename policies::normalise<
1037       Policy,
1038       policies::promote_float<false>,
1039       policies::promote_double<false>,
1040       policies::discrete_quantile<>,
1041       policies::assert_undefined<> >::type forwarding_policy;
1042
1043    if(a <= 0)
1044       return policies::raise_domain_error<result_type>(function, "The argument a to the incomplete beta function inverse must be greater than zero (got a=%1%).", a, pol);
1045    if(b <= 0)
1046       return policies::raise_domain_error<result_type>(function, "The argument b to the incomplete beta function inverse must be greater than zero (got b=%1%).", b, pol);
1047    if((q < 0) || (q > 1))
1048       return policies::raise_domain_error<result_type>(function, "Argument q outside the range [0,1] in the incomplete beta function inverse (got q=%1%).", q, pol);
1049
1050    value_type rx, ry;
1051
1052    rx = detail::ibeta_inv_imp(
1053          static_cast<value_type>(a),
1054          static_cast<value_type>(b),
1055          static_cast<value_type>(1 - q),
1056          static_cast<value_type>(q),
1057          forwarding_policy(), &ry);
1058
1059    if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function);
1060    return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function);
1061 }
1062
1063 template <class T1, class T2, class T3, class T4>
1064 inline typename tools::promote_args<T1, T2, T3, T4>::type
1065    ibetac_inv(T1 a, T2 b, T3 q, T4* py)
1066 {
1067    return ibetac_inv(a, b, q, py, policies::policy<>());
1068 }
1069
1070 template <class RT1, class RT2, class RT3>
1071 inline typename tools::promote_args<RT1, RT2, RT3>::type
1072    ibetac_inv(RT1 a, RT2 b, RT3 q)
1073 {
1074    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1075    return ibetac_inv(a, b, q, static_cast<result_type*>(nullptr), policies::policy<>());
1076 }
1077
1078 template <class RT1, class RT2, class RT3, class Policy>
1079 inline typename tools::promote_args<RT1, RT2, RT3>::type
1080    ibetac_inv(RT1 a, RT2 b, RT3 q, const Policy& pol)
1081 {
1082    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
1083    return ibetac_inv(a, b, q, static_cast<result_type*>(nullptr), pol);
1084 }
1085
1086 } // namespace math
1087 } // namespace boost
1088
1089 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
1090
1091
1092
1093