boost/math/special_functions/detail/ibeta_inverse.hpp

   1 //  Copyright John Maddock 2006.
   2 //  Copyright Paul A. Bristow 2007
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP
   8 #define BOOST_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP
   9
  10 #ifdef _MSC_VER
  11 #pragma once
  12 #endif
  13
  14 #include <boost/math/special_functions/beta.hpp>
  15 #include <boost/math/special_functions/erf.hpp>
  16 #include <boost/math/tools/roots.hpp>
  17 #include <boost/math/special_functions/detail/t_distribution_inv.hpp>
  18
  19 namespace boost{ namespace math{ namespace detail{
  20
  21 //
  22 // Helper object used by root finding
  23 // code to convert eta to x.
  24 //
  25 template <class T>
  26 struct temme_root_finder
  27 {
  28    temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {}
  29
  30    boost::math::tuple<T, T> operator()(T x)
  31    {
  32       BOOST_MATH_STD_USING // ADL of std names
  33
  34       T y = 1 - x;
  35       if(y == 0)
  36       {
  37          T big = tools::max_value<T>() / 4;
  38          return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big));
  39       }
  40       if(x == 0)
  41       {
  42          T big = tools::max_value<T>() / 4;
  43          return boost::math::make_tuple(static_cast<T>(-big), big);
  44       }
  45       T f = log(x) + a * log(y) + t;
  46       T f1 = (1 / x) - (a / (y));
  47       return boost::math::make_tuple(f, f1);
  48    }
  49 private:
  50    T t, a;
  51 };
  52 //
  53 // See:
  54 // "Asymptotic Inversion of the Incomplete Beta Function"
  55 // N.M. Temme
  56 // Journal of Computation and Applied Mathematics 41 (1992) 145-157.
  57 // Section 2.
  58 //
  59 template <class T, class Policy>
  60 T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)
  61 {
  62    BOOST_MATH_STD_USING // ADL of std names
  63
  64    const T r2 = sqrt(T(2));
  65    //
  66    // get the first approximation for eta from the inverse
  67    // error function (Eq: 2.9 and 2.10).
  68    //
  69    T eta0 = boost::math::erfc_inv(2 * z, pol);
  70    eta0 /= -sqrt(a / 2);
  71
  72    T terms[4] = { eta0 };
  73    T workspace[7];
  74    //
  75    // calculate powers:
  76    //
  77    T B = b - a;
  78    T B_2 = B * B;
  79    T B_3 = B_2 * B;
  80    //
  81    // Calculate correction terms:
  82    //
  83
  84    // See eq following 2.15:
  85    workspace[0] = -B * r2 / 2;
  86    workspace[1] = (1 - 2 * B) / 8;
  87    workspace[2] = -(B * r2 / 48);
  88    workspace[3] = T(-1) / 192;
  89    workspace[4] = -B * r2 / 3840;
  90    terms[1] = tools::evaluate_polynomial(workspace, eta0, 5);
  91    // Eq Following 2.17:
  92    workspace[0] = B * r2 * (3 * B - 2) / 12;
  93    workspace[1] = (20 * B_2 - 12 * B + 1) / 128;
  94    workspace[2] = B * r2 * (20 * B - 1) / 960;
  95    workspace[3] = (16 * B_2 + 30 * B - 15) / 4608;
  96    workspace[4] = B * r2 * (21 * B + 32) / 53760;
  97    workspace[5] = (-32 * B_2 + 63) / 368640;
  98    workspace[6] = -B * r2 * (120 * B + 17) / 25804480;
  99    terms[2] = tools::evaluate_polynomial(workspace, eta0, 7);
 100    // Eq Following 2.17:
 101    workspace[0] = B * r2 * (-75 * B_2 + 80 * B - 16) / 480;
 102    workspace[1] = (-1080 * B_3 + 868 * B_2 - 90 * B - 45) / 9216;
 103    workspace[2] = B * r2 * (-1190 * B_2 + 84 * B + 373) / 53760;
 104    workspace[3] = (-2240 * B_3 - 2508 * B_2 + 2100 * B - 165) / 368640;
 105    terms[3] = tools::evaluate_polynomial(workspace, eta0, 4);
 106    //
 107    // Bring them together to get a final estimate for eta:
 108    //
 109    T eta = tools::evaluate_polynomial(terms, T(1/a), 4);
 110    //
 111    // now we need to convert eta to x, by solving the appropriate
 112    // quadratic equation:
 113    //
 114    T eta_2 = eta * eta;
 115    T c = -exp(-eta_2 / 2);
 116    T x;
 117    if(eta_2 == 0)
 118       x = 0.5;
 119    else
 120       x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;
 121
 122    BOOST_ASSERT(x >= 0);
 123    BOOST_ASSERT(x <= 1);
 124    BOOST_ASSERT(eta * (x - 0.5) >= 0);
 125 #ifdef BOOST_INSTRUMENT
 126    std::cout << "Estimating x with Temme method 1: " << x << std::endl;
 127 #endif
 128    return x;
 129 }
 130 //
 131 // See:
 132 // "Asymptotic Inversion of the Incomplete Beta Function"
 133 // N.M. Temme
 134 // Journal of Computation and Applied Mathematics 41 (1992) 145-157.
 135 // Section 3.
 136 //
 137 template <class T, class Policy>
 138 T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r, T theta, const Policy& pol)
 139 {
 140    BOOST_MATH_STD_USING // ADL of std names
 141
 142    //
 143    // Get first estimate for eta, see Eq 3.9 and 3.10,
 144    // but note there is a typo in Eq 3.10:
 145    //
 146    T eta0 = boost::math::erfc_inv(2 * z, pol);
 147    eta0 /= -sqrt(r / 2);
 148
 149    T s = sin(theta);
 150    T c = cos(theta);
 151    //
 152    // Now we need to purturb eta0 to get eta, which we do by
 153    // evaluating the polynomial in 1/r at the bottom of page 151,
 154    // to do this we first need the error terms e1, e2 e3
 155    // which we'll fill into the array "terms".  Since these
 156    // terms are themselves polynomials, we'll need another
 157    // array "workspace" to calculate those...
 158    //
 159    T terms[4] = { eta0 };
 160    T workspace[6];
 161    //
 162    // some powers of sin(theta)cos(theta) that we'll need later:
 163    //
 164    T sc = s * c;
 165    T sc_2 = sc * sc;
 166    T sc_3 = sc_2 * sc;
 167    T sc_4 = sc_2 * sc_2;
 168    T sc_5 = sc_2 * sc_3;
 169    T sc_6 = sc_3 * sc_3;
 170    T sc_7 = sc_4 * sc_3;
 171    //
 172    // Calculate e1 and put it in terms[1], see the middle of page 151:
 173    //
 174    workspace[0] = (2 * s * s - 1) / (3 * s * c);
 175    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co1[] = { -1, -5, 5 };
 176    workspace[1] = -tools::evaluate_even_polynomial(co1, s, 3) / (36 * sc_2);
 177    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co2[] = { 1, 21, -69, 46 };
 178    workspace[2] = tools::evaluate_even_polynomial(co2, s, 4) / (1620 * sc_3);
 179    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co3[] = { 7, -2, 33, -62, 31 };
 180    workspace[3] = -tools::evaluate_even_polynomial(co3, s, 5) / (6480 * sc_4);
 181    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co4[] = { 25, -52, -17, 88, -115, 46 };
 182    workspace[4] = tools::evaluate_even_polynomial(co4, s, 6) / (90720 * sc_5);
 183    terms[1] = tools::evaluate_polynomial(workspace, eta0, 5);
 184    //
 185    // Now evaluate e2 and put it in terms[2]:
 186    //
 187    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co5[] = { 7, 12, -78, 52 };
 188    workspace[0] = -tools::evaluate_even_polynomial(co5, s, 4) / (405 * sc_3);
 189    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co6[] = { -7, 2, 183, -370, 185 };
 190    workspace[1] = tools::evaluate_even_polynomial(co6, s, 5) / (2592 * sc_4);
 191    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co7[] = { -533, 776, -1835, 10240, -13525, 5410 };
 192    workspace[2] = -tools::evaluate_even_polynomial(co7, s, 6) / (204120 * sc_5);
 193    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co8[] = { -1579, 3747, -3372, -15821, 45588, -45213, 15071 };
 194    workspace[3] = -tools::evaluate_even_polynomial(co8, s, 7) / (2099520 * sc_6);
 195    terms[2] = tools::evaluate_polynomial(workspace, eta0, 4);
 196    //
 197    // And e3, and put it in terms[3]:
 198    //
 199    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co9[] = {449, -1259, -769, 6686, -9260, 3704 };
 200    workspace[0] = tools::evaluate_even_polynomial(co9, s, 6) / (102060 * sc_5);
 201    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co10[] = { 63149, -151557, 140052, -727469, 2239932, -2251437, 750479 };
 202    workspace[1] = -tools::evaluate_even_polynomial(co10, s, 7) / (20995200 * sc_6);
 203    static const BOOST_MATH_INT_TABLE_TYPE(T, int) co11[] = { 29233, -78755, 105222, 146879, -1602610, 3195183, -2554139, 729754 };
 204    workspace[2] = tools::evaluate_even_polynomial(co11, s, 8) / (36741600 * sc_7);
 205    terms[3] = tools::evaluate_polynomial(workspace, eta0, 3);
 206    //
 207    // Bring the correction terms together to evaluate eta,
 208    // this is the last equation on page 151:
 209    //
 210    T eta = tools::evaluate_polynomial(terms, T(1/r), 4);
 211    //
 212    // Now that we have eta we need to back solve for x,
 213    // we seek the value of x that gives eta in Eq 3.2.
 214    // The two methods used are described in section 5.
 215    //
 216    // Begin by defining a few variables we'll need later:
 217    //
 218    T x;
 219    T s_2 = s * s;
 220    T c_2 = c * c;
 221    T alpha = c / s;
 222    alpha *= alpha;
 223    T lu = (-(eta * eta) / (2 * s_2) + log(s_2) + c_2 * log(c_2) / s_2);
 224    //
 225    // Temme doesn't specify what value to switch on here,
 226    // but this seems to work pretty well:
 227    //
 228    if(fabs(eta) < 0.7)
 229    {
 230       //
 231       // Small eta use the expansion Temme gives in the second equation
 232       // of section 5, it's a polynomial in eta:
 233       //
 234       workspace[0] = s * s;
 235       workspace[1] = s * c;
 236       workspace[2] = (1 - 2 * workspace[0]) / 3;
 237       static const BOOST_MATH_INT_TABLE_TYPE(T, int) co12[] = { 1, -13, 13 };
 238       workspace[3] = tools::evaluate_polynomial(co12, workspace[0], 3) / (36 * s * c);
 239       static const BOOST_MATH_INT_TABLE_TYPE(T, int) co13[] = { 1, 21, -69, 46 };
 240       workspace[4] = tools::evaluate_polynomial(co13, workspace[0], 4) / (270 * workspace[0] * c * c);
 241       x = tools::evaluate_polynomial(workspace, eta, 5);
 242 #ifdef BOOST_INSTRUMENT
 243       std::cout << "Estimating x with Temme method 2 (small eta): " << x << std::endl;
 244 #endif
 245    }
 246    else
 247    {
 248       //
 249       // If eta is large we need to solve Eq 3.2 more directly,
 250       // begin by getting an initial approximation for x from
 251       // the last equation on page 155, this is a polynomial in u:
 252       //
 253       T u = exp(lu);
 254       workspace[0] = u;
 255       workspace[1] = alpha;
 256       workspace[2] = 0;
 257       workspace[3] = 3 * alpha * (3 * alpha + 1) / 6;
 258       workspace[4] = 4 * alpha * (4 * alpha + 1) * (4 * alpha + 2) / 24;
 259       workspace[5] = 5 * alpha * (5 * alpha + 1) * (5 * alpha + 2) * (5 * alpha + 3) / 120;
 260       x = tools::evaluate_polynomial(workspace, u, 6);
 261       //
 262       // At this point we may or may not have the right answer, Eq-3.2 has
 263       // two solutions for x for any given eta, however the mapping in 3.2
 264       // is 1:1 with the sign of eta and x-sin^2(theta) being the same.
 265       // So we can check if we have the right root of 3.2, and if not
 266       // switch x for 1-x.  This transformation is motivated by the fact
 267       // that the distribution is *almost* symetric so 1-x will be in the right
 268       // ball park for the solution:
 269       //
 270       if((x - s_2) * eta < 0)
 271          x = 1 - x;
 272 #ifdef BOOST_INSTRUMENT
 273       std::cout << "Estimating x with Temme method 2 (large eta): " << x << std::endl;
 274 #endif
 275    }
 276    //
 277    // The final step is a few Newton-Raphson iterations to
 278    // clean up our approximation for x, this is pretty cheap
 279    // in general, and very cheap compared to an incomplete beta
 280    // evaluation.  The limits set on x come from the observation
 281    // that the sign of eta and x-sin^2(theta) are the same.
 282    //
 283    T lower, upper;
 284    if(eta < 0)
 285    {
 286       lower = 0;
 287       upper = s_2;
 288    }
 289    else
 290    {
 291       lower = s_2;
 292       upper = 1;
 293    }
 294    //
 295    // If our initial approximation is out of bounds then bisect:
 296    //
 297    if((x < lower) || (x > upper))
 298       x = (lower+upper) / 2;
 299    //
 300    // And iterate:
 301    //
 302    x = tools::newton_raphson_iterate(
 303       temme_root_finder<T>(-lu, alpha), x, lower, upper, policies::digits<T, Policy>() / 2);
 304
 305    return x;
 306 }
 307 //
 308 // See:
 309 // "Asymptotic Inversion of the Incomplete Beta Function"
 310 // N.M. Temme
 311 // Journal of Computation and Applied Mathematics 41 (1992) 145-157.
 312 // Section 4.
 313 //
 314 template <class T, class Policy>
 315 T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
 316 {
 317    BOOST_MATH_STD_USING // ADL of std names
 318
 319    //
 320    // Begin by getting an initial approximation for the quantity
 321    // eta from the dominant part of the incomplete beta:
 322    //
 323    T eta0;
 324    if(p < q)
 325       eta0 = boost::math::gamma_q_inv(b, p, pol);
 326    else
 327       eta0 = boost::math::gamma_p_inv(b, q, pol);
 328    eta0 /= a;
 329    //
 330    // Define the variables and powers we'll need later on:
 331    //
 332    T mu = b / a;
 333    T w = sqrt(1 + mu);
 334    T w_2 = w * w;
 335    T w_3 = w_2 * w;
 336    T w_4 = w_2 * w_2;
 337    T w_5 = w_3 * w_2;
 338    T w_6 = w_3 * w_3;
 339    T w_7 = w_4 * w_3;
 340    T w_8 = w_4 * w_4;
 341    T w_9 = w_5 * w_4;
 342    T w_10 = w_5 * w_5;
 343    T d = eta0 - mu;
 344    T d_2 = d * d;
 345    T d_3 = d_2 * d;
 346    T d_4 = d_2 * d_2;
 347    T w1 = w + 1;
 348    T w1_2 = w1 * w1;
 349    T w1_3 = w1 * w1_2;
 350    T w1_4 = w1_2 * w1_2;
 351    //
 352    // Now we need to compute the purturbation error terms that
 353    // convert eta0 to eta, these are all polynomials of polynomials.
 354    // Probably these should be re-written to use tabulated data
 355    // (see examples above), but it's less of a win in this case as we
 356    // need to calculate the individual powers for the denominator terms
 357    // anyway, so we might as well use them for the numerator-polynomials
 358    // as well....
 359    //
 360    // Refer to p154-p155 for the details of these expansions:
 361    //
 362    T e1 = (w + 2) * (w - 1) / (3 * w);
 363    e1 += (w_3 + 9 * w_2 + 21 * w + 5) * d / (36 * w_2 * w1);
 364    e1 -= (w_4 - 13 * w_3 + 69 * w_2 + 167 * w + 46) * d_2 / (1620 * w1_2 * w_3);
 365    e1 -= (7 * w_5 + 21 * w_4 + 70 * w_3 + 26 * w_2 - 93 * w - 31) * d_3 / (6480 * w1_3 * w_4);
 366    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);
 367
 368    T e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1) / (1620 * w1 * w_3);
 369    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);
 370    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);
 371    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 / (14696640 * w1_4 * w_6);
 372
 373    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);
 374    e3 -= (442043 * w_9 + 2054169 * w_8 + 3803094 * w_7 + 3470754 * w_6 + 2141568 * w_5 - 2393568 * w_4 - 19904934 * w_3 - 34714674 * w_2 - 23128299 * w - 5253353) * d / (146966400 * w_6 * w1_3);
 375    e3 -= (116932 * w_10 + 819281 * w_9 + 2378172 * w_8 + 4341330 * w_7 + 6806004 * w_6 + 10622748 * w_5 + 18739500 * w_4 + 30651894 * w_3 + 30869976 * w_2 + 15431867 * w + 2919016) * d_2 / (146966400 * w1_4 * w_7);
 376    //
 377    // Combine eta0 and the error terms to compute eta (Second eqaution p155):
 378    //
 379    T eta = eta0 + e1 / a + e2 / (a * a) + e3 / (a * a * a);
 380    //
 381    // Now we need to solve Eq 4.2 to obtain x.  For any given value of
 382    // eta there are two solutions to this equation, and since the distribtion
 383    // may be very skewed, these are not related by x ~ 1-x we used when
 384    // implementing section 3 above.  However we know that:
 385    //
 386    //  cross < x <= 1       ; iff eta < mu
 387    //          x == cross   ; iff eta == mu
 388    //     0 <= x < cross    ; iff eta > mu
 389    //
 390    // Where cross == 1 / (1 + mu)
 391    // Many thanks to Prof Temme for clarifying this point.
 392    //
 393    // Therefore we'll just jump straight into Newton iterations
 394    // to solve Eq 4.2 using these bounds, and simple bisection
 395    // as the first guess, in practice this converges pretty quickly
 396    // and we only need a few digits correct anyway:
 397    //
 398    if(eta <= 0)
 399       eta = tools::min_value<T>();
 400    T u = eta - mu * log(eta) + (1 + mu) * log(1 + mu) - mu;
 401    T cross = 1 / (1 + mu);
 402    T lower = eta < mu ? cross : 0;
 403    T upper = eta < mu ? 1 : cross;
 404    T x = (lower + upper) / 2;
 405    x = tools::newton_raphson_iterate(
 406       temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2);
 407 #ifdef BOOST_INSTRUMENT
 408    std::cout << "Estimating x with Temme method 3: " << x << std::endl;
 409 #endif
 410    return x;
 411 }
 412
 413 template <class T, class Policy>
 414 struct ibeta_roots
 415 {
 416    ibeta_roots(T _a, T _b, T t, bool inv = false)
 417       : a(_a), b(_b), target(t), invert(inv) {}
 418
 419    boost::math::tuple<T, T, T> operator()(T x)
 420    {
 421       BOOST_MATH_STD_USING // ADL of std names
 422
 423       BOOST_FPU_EXCEPTION_GUARD
 424
 425       T f1;
 426       T y = 1 - x;
 427       T f = ibeta_imp(a, b, x, Policy(), invert, true, &f1) - target;
 428       if(invert)
 429          f1 = -f1;
 430       if(y == 0)
 431          y = tools::min_value<T>() * 64;
 432       if(x == 0)
 433          x = tools::min_value<T>() * 64;
 434
 435       T f2 = f1 * (-y * a + (b - 2) * x + 1);
 436       if(fabs(f2) < y * x * tools::max_value<T>())
 437          f2 /= (y * x);
 438       if(invert)
 439          f2 = -f2;
 440
 441       // make sure we don't have a zero derivative:
 442       if(f1 == 0)
 443          f1 = (invert ? -1 : 1) * tools::min_value<T>() * 64;
 444
 445       return boost::math::make_tuple(f, f1, f2);
 446    }
 447 private:
 448    T a, b, target;
 449    bool invert;
 450 };
 451
 452 template <class T, class Policy>
 453 T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
 454 {
 455    BOOST_MATH_STD_USING  // For ADL of math functions.
 456
 457    //
 458    // The flag invert is set to true if we swap a for b and p for q,
 459    // in which case the result has to be subtracted from 1:
 460    //
 461    bool invert = false;
 462    //
 463    // Handle trivial cases first:
 464    //
 465    if(q == 0)
 466    {
 467       if(py) *py = 0;
 468       return 1;
 469    }
 470    else if(p == 0)
 471    {
 472       if(py) *py = 1;
 473       return 0;
 474    }
 475    else if(a == 1)
 476    {
 477       if(b == 1)
 478       {
 479          if(py) *py = 1 - p;
 480          return p;
 481       }
 482       // Change things around so we can handle as b == 1 special case below:
 483       std::swap(a, b);
 484       std::swap(p, q);
 485       invert = true;
 486    }
 487    //
 488    // Depending upon which approximation method we use, we may end up
 489    // calculating either x or y initially (where y = 1-x):
 490    //
 491    T x = 0; // Set to a safe zero to avoid a
 492    // MSVC 2005 warning C4701: potentially uninitialized local variable 'x' used
 493    // But code inspection appears to ensure that x IS assigned whatever the code path.
 494    T y;
 495
 496    // For some of the methods we can put tighter bounds
 497    // on the result than simply [0,1]:
 498    //
 499    T lower = 0;
 500    T upper = 1;
 501    //
 502    // Student's T with b = 0.5 gets handled as a special case, swap
 503    // around if the arguments are in the "wrong" order:
 504    //
 505    if(a == 0.5f)
 506    {
 507       if(b == 0.5f)
 508       {
 509          x = sin(p * constants::half_pi<T>());
 510          x *= x;
 511          if(py)
 512          {
 513             *py = sin(q * constants::half_pi<T>());
 514             *py *= *py;
 515          }
 516          return x;
 517       }
 518       else if(b > 0.5f)
 519       {
 520          std::swap(a, b);
 521          std::swap(p, q);
 522          invert = !invert;
 523       }
 524    }
 525    //
 526    // Select calculation method for the initial estimate:
 527    //
 528    if((b == 0.5f) && (a >= 0.5f) && (p != 1))
 529    {
 530       //
 531       // We have a Student's T distribution:
 532       x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol);
 533    }
 534    else if(b == 1)
 535    {
 536       if(p < q)
 537       {
 538          if(a > 1)
 539          {
 540             x = pow(p, 1 / a);
 541             y = -expm1(log(p) / a);
 542          }
 543          else
 544          {
 545             x = pow(p, 1 / a);
 546             y = 1 - x;
 547          }
 548       }
 549       else
 550       {
 551          x = exp(log1p(-q) / a);
 552          y = -expm1(log1p(-q) / a);
 553       }
 554       if(invert)
 555          std::swap(x, y);
 556       if(py)
 557          *py = y;
 558       return x;
 559    }
 560    else if(a + b > 5)
 561    {
 562       //
 563       // When a+b is large then we can use one of Prof Temme's
 564       // asymptotic expansions, begin by swapping things around
 565       // so that p < 0.5, we do this to avoid cancellations errors
 566       // when p is large.
 567       //
 568       if(p > 0.5)
 569       {
 570          std::swap(a, b);
 571          std::swap(p, q);
 572          invert = !invert;
 573       }
 574       T minv = (std::min)(a, b);
 575       T maxv = (std::max)(a, b);
 576       if((sqrt(minv) > (maxv - minv)) && (minv > 5))
 577       {
 578          //
 579          // When a and b differ by a small amount
 580          // the curve is quite symmetrical and we can use an error
 581          // function to approximate the inverse. This is the cheapest
 582          // of the three Temme expantions, and the calculated value
 583          // for x will never be much larger than p, so we don't have
 584          // to worry about cancellation as long as p is small.
 585          //
 586          x = temme_method_1_ibeta_inverse(a, b, p, pol);
 587          y = 1 - x;
 588       }
 589       else
 590       {
 591          T r = a + b;
 592          T theta = asin(sqrt(a / r));
 593          T lambda = minv / r;
 594          if((lambda >= 0.2) && (lambda <= 0.8) && (r >= 10))
 595          {
 596             //
 597             // The second error function case is the next cheapest
 598             // to use, it brakes down when the result is likely to be
 599             // very small, if a+b is also small, but we can use a
 600             // cheaper expansion there in any case.  As before x won't
 601             // be much larger than p, so as long as p is small we should
 602             // be free of cancellation error.
 603             //
 604             T ppa = pow(p, 1/a);
 605             if((ppa < 0.0025) && (a + b < 200))
 606             {
 607                x = ppa * pow(a * boost::math::beta(a, b, pol), 1/a);
 608             }
 609             else
 610                x = temme_method_2_ibeta_inverse(a, b, p, r, theta, pol);
 611             y = 1 - x;
 612          }
 613          else
 614          {
 615             //
 616             // If we get here then a and b are very different in magnitude
 617             // and we need to use the third of Temme's methods which
 618             // involves inverting the incomplete gamma.  This is much more
 619             // expensive than the other methods.  We also can only use this
 620             // method when a > b, which can lead to cancellation errors
 621             // if we really want y (as we will when x is close to 1), so
 622             // a different expansion is used in that case.
 623             //
 624             if(a < b)
 625             {
 626                std::swap(a, b);
 627                std::swap(p, q);
 628                invert = !invert;
 629             }
 630             //
 631             // Try and compute the easy way first:
 632             //
 633             T bet = 0;
 634             if(b < 2)
 635                bet = boost::math::beta(a, b, pol);
 636             if(bet != 0)
 637             {
 638                y = pow(b * q * bet, 1/b);
 639                x = 1 - y;
 640             }
 641             else
 642                y = 1;
 643             if(y > 1e-5)
 644             {
 645                x = temme_method_3_ibeta_inverse(a, b, p, q, pol);
 646                y = 1 - x;
 647             }
 648          }
 649       }
 650    }
 651    else if((a < 1) && (b < 1))
 652    {
 653       //
 654       // Both a and b less than 1,
 655       // there is a point of inflection at xs:
 656       //
 657       T xs = (1 - a) / (2 - a - b);
 658       //
 659       // Now we need to ensure that we start our iteration from the
 660       // right side of the inflection point:
 661       //
 662       T fs = boost::math::ibeta(a, b, xs, pol) - p;
 663       if(fabs(fs) / p < tools::epsilon<T>() * 3)
 664       {
 665          // The result is at the point of inflection, best just return it:
 666          *py = invert ? xs : 1 - xs;
 667          return invert ? 1-xs : xs;
 668       }
 669       if(fs < 0)
 670       {
 671          std::swap(a, b);
 672          std::swap(p, q);
 673          invert = !invert;
 674          xs = 1 - xs;
 675       }
 676       T xg = pow(a * p * boost::math::beta(a, b, pol), 1/a);
 677       x = xg / (1 + xg);
 678       y = 1 / (1 + xg);
 679       //
 680       // And finally we know that our result is below the inflection
 681       // point, so set an upper limit on our search:
 682       //
 683       if(x > xs)
 684          x = xs;
 685       upper = xs;
 686    }
 687    else if((a > 1) && (b > 1))
 688    {
 689       //
 690       // Small a and b, both greater than 1,
 691       // there is a point of inflection at xs,
 692       // and it's complement is xs2, we must always
 693       // start our iteration from the right side of the
 694       // point of inflection.
 695       //
 696       T xs = (a - 1) / (a + b - 2);
 697       T xs2 = (b - 1) / (a + b - 2);
 698       T ps = boost::math::ibeta(a, b, xs, pol) - p;
 699
 700       if(ps < 0)
 701       {
 702          std::swap(a, b);
 703          std::swap(p, q);
 704          std::swap(xs, xs2);
 705          invert = !invert;
 706       }
 707       //
 708       // Estimate x and y, using expm1 to get a good estimate
 709       // for y when it's very small:
 710       //
 711       T lx = log(p * a * boost::math::beta(a, b, pol)) / a;
 712       x = exp(lx);
 713       y = x < 0.9 ? T(1 - x) : (T)(-boost::math::expm1(lx, pol));
 714
 715       if((b < a) && (x < 0.2))
 716       {
 717          //
 718          // Under a limited range of circumstances we can improve
 719          // our estimate for x, frankly it's clear if this has much effect!
 720          //
 721          T ap1 = a - 1;
 722          T bm1 = b - 1;
 723          T a_2 = a * a;
 724          T a_3 = a * a_2;
 725          T b_2 = b * b;
 726          T terms[5] = { 0, 1 };
 727          terms[2] = bm1 / ap1;
 728          ap1 *= ap1;
 729          terms[3] = bm1 * (3 * a * b + 5 * b + a_2 - a - 4) / (2 * (a + 2) * ap1);
 730          ap1 *= (a + 1);
 731          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)
 732                     / (3 * (a + 3) * (a + 2) * ap1);
 733          x = tools::evaluate_polynomial(terms, x, 5);
 734       }
 735       //
 736       // And finally we know that our result is below the inflection
 737       // point, so set an upper limit on our search:
 738       //
 739       if(x > xs)
 740          x = xs;
 741       upper = xs;
 742    }
 743    else /*if((a <= 1) != (b <= 1))*/
 744    {
 745       //
 746       // If all else fails we get here, only one of a and b
 747       // is above 1, and a+b is small.  Start by swapping
 748       // things around so that we have a concave curve with b > a
 749       // and no points of inflection in [0,1].  As long as we expect
 750       // x to be small then we can use the simple (and cheap) power
 751       // term to estimate x, but when we expect x to be large then
 752       // this greatly underestimates x and leaves us trying to
 753       // iterate "round the corner" which may take almost forever...
 754       //
 755       // We could use Temme's inverse gamma function case in that case,
 756       // this works really rather well (albeit expensively) even though
 757       // strictly speaking we're outside it's defined range.
 758       //
 759       // However it's expensive to compute, and an alternative approach
 760       // which models the curve as a distorted quarter circle is much
 761       // cheaper to compute, and still keeps the number of iterations
 762       // required down to a reasonable level.  With thanks to Prof Temme
 763       // for this suggestion.
 764       //
 765       if(b < a)
 766       {
 767          std::swap(a, b);
 768          std::swap(p, q);
 769          invert = !invert;
 770       }
 771       if(pow(p, 1/a) < 0.5)
 772       {
 773          x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
 774          if(x == 0)
 775             x = boost::math::tools::min_value<T>();
 776          y = 1 - x;
 777       }
 778       else /*if(pow(q, 1/b) < 0.1)*/
 779       {
 780          // model a distorted quarter circle:
 781          y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
 782          if(y == 0)
 783             y = boost::math::tools::min_value<T>();
 784          x = 1 - y;
 785       }
 786    }
 787
 788    //
 789    // Now we have a guess for x (and for y) we can set things up for
 790    // iteration.  If x > 0.5 it pays to swap things round:
 791    //
 792    if(x > 0.5)
 793    {
 794       std::swap(a, b);
 795       std::swap(p, q);
 796       std::swap(x, y);
 797       invert = !invert;
 798       T l = 1 - upper;
 799       T u = 1 - lower;
 800       lower = l;
 801       upper = u;
 802    }
 803    //
 804    // lower bound for our search:
 805    //
 806    // We're not interested in denormalised answers as these tend to
 807    // these tend to take up lots of iterations, given that we can't get
 808    // accurate derivatives in this area (they tend to be infinite).
 809    //
 810    if(lower == 0)
 811    {
 812       if(invert && (py == 0))
 813       {
 814          //
 815          // We're not interested in answers smaller than machine epsilon:
 816          //
 817          lower = boost::math::tools::epsilon<T>();
 818          if(x < lower)
 819             x = lower;
 820       }
 821       else
 822          lower = boost::math::tools::min_value<T>();
 823       if(x < lower)
 824          x = lower;
 825    }
 826    //
 827    // Figure out how many digits to iterate towards:
 828    //
 829    int digits = boost::math::policies::digits<T, Policy>() / 2;
 830    if((x < 1e-50) && ((a < 1) || (b < 1)))
 831    {
 832       //
 833       // If we're in a region where the first derivative is very
 834       // large, then we have to take care that the root-finder
 835       // doesn't terminate prematurely.  We'll bump the precision
 836       // up to avoid this, but we have to take care not to set the
 837       // precision too high or the last few iterations will just
 838       // thrash around and convergence may be slow in this case.
 839       // Try 3/4 of machine epsilon:
 840       //
 841       digits *= 3;
 842       digits /= 2;
 843    }
 844    //
 845    // Now iterate, we can use either p or q as the target here
 846    // depending on which is smaller:
 847    //
 848    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 849    x = boost::math::tools::halley_iterate(
 850       boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
 851    policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
 852    //
 853    // We don't really want these asserts here, but they are useful for sanity
 854    // checking that we have the limits right, uncomment if you suspect bugs *only*.
 855    //
 856    //BOOST_ASSERT(x != upper);
 857    //BOOST_ASSERT((x != lower) || (x == boost::math::tools::min_value<T>()) || (x == boost::math::tools::epsilon<T>()));
 858    //
 859    // Tidy up, if we "lower" was too high then zero is the best answer we have:
 860    //
 861    if(x == lower)
 862       x = 0;
 863    if(py)
 864       *py = invert ? x : 1 - x;
 865    return invert ? 1-x : x;
 866 }
 867
 868 } // namespace detail
 869
 870 template <class T1, class T2, class T3, class T4, class Policy>
 871 inline typename tools::promote_args<T1, T2, T3, T4>::type
 872    ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy& pol)
 873 {
 874    static const char* function = "boost::math::ibeta_inv<%1%>(%1%,%1%,%1%)";
 875    BOOST_FPU_EXCEPTION_GUARD
 876    typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type;
 877    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 878    typedef typename policies::normalise<
 879       Policy,
 880       policies::promote_float<false>,
 881       policies::promote_double<false>,
 882       policies::discrete_quantile<>,
 883       policies::assert_undefined<> >::type forwarding_policy;
 884
 885    if(a <= 0)
 886       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);
 887    if(b <= 0)
 888       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);
 889    if((p < 0) || (p > 1))
 890       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);
 891
 892    value_type rx, ry;
 893
 894    rx = detail::ibeta_inv_imp(
 895          static_cast<value_type>(a),
 896          static_cast<value_type>(b),
 897          static_cast<value_type>(p),
 898          static_cast<value_type>(1 - p),
 899          forwarding_policy(), &ry);
 900
 901    if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function);
 902    return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function);
 903 }
 904
 905 template <class T1, class T2, class T3, class T4>
 906 inline typename tools::promote_args<T1, T2, T3, T4>::type
 907    ibeta_inv(T1 a, T2 b, T3 p, T4* py)
 908 {
 909    return ibeta_inv(a, b, p, py, policies::policy<>());
 910 }
 911
 912 template <class T1, class T2, class T3>
 913 inline typename tools::promote_args<T1, T2, T3>::type
 914    ibeta_inv(T1 a, T2 b, T3 p)
 915 {
 916    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
 917    return ibeta_inv(a, b, p, static_cast<result_type*>(0), policies::policy<>());
 918 }
 919
 920 template <class T1, class T2, class T3, class Policy>
 921 inline typename tools::promote_args<T1, T2, T3>::type
 922    ibeta_inv(T1 a, T2 b, T3 p, const Policy& pol)
 923 {
 924    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
 925    return ibeta_inv(a, b, p, static_cast<result_type*>(0), pol);
 926 }
 927
 928 template <class T1, class T2, class T3, class T4, class Policy>
 929 inline typename tools::promote_args<T1, T2, T3, T4>::type
 930    ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy& pol)
 931 {
 932    static const char* function = "boost::math::ibetac_inv<%1%>(%1%,%1%,%1%)";
 933    BOOST_FPU_EXCEPTION_GUARD
 934    typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type;
 935    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 936    typedef typename policies::normalise<
 937       Policy,
 938       policies::promote_float<false>,
 939       policies::promote_double<false>,
 940       policies::discrete_quantile<>,
 941       policies::assert_undefined<> >::type forwarding_policy;
 942
 943    if(a <= 0)
 944       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);
 945    if(b <= 0)
 946       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);
 947    if((q < 0) || (q > 1))
 948       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);
 949
 950    value_type rx, ry;
 951
 952    rx = detail::ibeta_inv_imp(
 953          static_cast<value_type>(a),
 954          static_cast<value_type>(b),
 955          static_cast<value_type>(1 - q),
 956          static_cast<value_type>(q),
 957          forwarding_policy(), &ry);
 958
 959    if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function);
 960    return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function);
 961 }
 962
 963 template <class T1, class T2, class T3, class T4>
 964 inline typename tools::promote_args<T1, T2, T3, T4>::type
 965    ibetac_inv(T1 a, T2 b, T3 q, T4* py)
 966 {
 967    return ibetac_inv(a, b, q, py, policies::policy<>());
 968 }
 969
 970 template <class RT1, class RT2, class RT3>
 971 inline typename tools::promote_args<RT1, RT2, RT3>::type
 972    ibetac_inv(RT1 a, RT2 b, RT3 q)
 973 {
 974    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 975    return ibetac_inv(a, b, q, static_cast<result_type*>(0), policies::policy<>());
 976 }
 977
 978 template <class RT1, class RT2, class RT3, class Policy>
 979 inline typename tools::promote_args<RT1, RT2, RT3>::type
 980    ibetac_inv(RT1 a, RT2 b, RT3 q, const Policy& pol)
 981 {
 982    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 983    return ibetac_inv(a, b, q, static_cast<result_type*>(0), pol);
 984 }
 985
 986 } // namespace math
 987 } // namespace boost
 988
 989 #endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP
 990
 991
 992
 993