1123 lines · plain
1// Copyright John Maddock 2006.2// Copyright Paul A. Bristow 20073// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP8#define BOOST_MATH_SPECIAL_FUNCTIONS_IBETA_INVERSE_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/precision.hpp>16#include <boost/math/tools/roots.hpp>17#include <boost/math/tools/tuple.hpp>18#include <boost/math/special_functions/beta.hpp>19#include <boost/math/special_functions/erf.hpp>20#include <boost/math/special_functions/detail/t_distribution_inv.hpp>21#include <boost/math/special_functions/fpclassify.hpp>22 23namespace boost{ namespace math{ namespace detail{24 25//26// Helper object used by root finding27// code to convert eta to x.28//29template <class T>30struct temme_root_finder31{32 BOOST_MATH_GPU_ENABLED temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {33 BOOST_MATH_ASSERT(34 math::tools::epsilon<T>() <= a && !(boost::math::isinf)(a));35 }36 37 BOOST_MATH_GPU_ENABLED boost::math::tuple<T, T> operator()(T x)38 {39 BOOST_MATH_STD_USING // ADL of std names40 41 T y = 1 - x;42 T f = log(x) + a * log(y) + t;43 T f1 = (1 / x) - (a / (y));44 return boost::math::make_tuple(f, f1);45 }46private:47 T t, a;48};49//50// See:51// "Asymptotic Inversion of the Incomplete Beta Function"52// N.M. Temme53// Journal of Computation and Applied Mathematics 41 (1992) 145-157.54// Section 2.55//56template <class T, class Policy>57BOOST_MATH_GPU_ENABLED T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)58{59 BOOST_MATH_STD_USING // ADL of std names60 61 const T r2 = sqrt(T(2));62 //63 // get the first approximation for eta from the inverse64 // error function (Eq: 2.9 and 2.10).65 //66 T eta0 = boost::math::erfc_inv(2 * z, pol);67 eta0 /= -sqrt(a / 2);68 69 T terms[4] = { eta0 };70 T workspace[7];71 //72 // calculate powers:73 //74 T B = b - a;75 T B_2 = B * B;76 T B_3 = B_2 * B;77 //78 // Calculate correction terms:79 //80 81 // See eq following 2.15:82 workspace[0] = -B * r2 / 2;83 workspace[1] = (1 - 2 * B) / 8;84 workspace[2] = -(B * r2 / 48);85 workspace[3] = T(-1) / 192;86 workspace[4] = -B * r2 / 3840;87 terms[1] = tools::evaluate_polynomial(workspace, eta0, 5);88 // Eq Following 2.17:89 workspace[0] = B * r2 * (3 * B - 2) / 12;90 workspace[1] = (20 * B_2 - 12 * B + 1) / 128;91 workspace[2] = B * r2 * (20 * B - 1) / 960;92 workspace[3] = (16 * B_2 + 30 * B - 15) / 4608;93 workspace[4] = B * r2 * (21 * B + 32) / 53760;94 workspace[5] = (-32 * B_2 + 63) / 368640;95 workspace[6] = -B * r2 * (120 * B + 17) / 25804480;96 terms[2] = tools::evaluate_polynomial(workspace, eta0, 7);97 // Eq Following 2.17:98 workspace[0] = B * r2 * (-75 * B_2 + 80 * B - 16) / 480;99 workspace[1] = (-1080 * B_3 + 868 * B_2 - 90 * B - 45) / 9216;100 workspace[2] = B * r2 * (-1190 * B_2 + 84 * B + 373) / 53760;101 workspace[3] = (-2240 * B_3 - 2508 * B_2 + 2100 * B - 165) / 368640;102 terms[3] = tools::evaluate_polynomial(workspace, eta0, 4);103 //104 // Bring them together to get a final estimate for eta:105 //106 T eta = tools::evaluate_polynomial(terms, T(1/a), 4);107 //108 // now we need to convert eta to x, by solving the appropriate109 // quadratic equation:110 //111 T eta_2 = eta * eta;112 T c = -exp(-eta_2 / 2);113 T x;114 if(eta_2 == 0)115 x = static_cast<T>(0.5f);116 else117 x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;118 //119 // These are post-conditions of the method, but the addition above120 // may result in us being out by 1ulp either side of the boundary,121 // so just check that we're in bounds and adjust as needed.122 // See https://github.com/boostorg/math/issues/961123 //124 if (x < 0)125 x = 0;126 else if (x > 1)127 x = 1;128 129 BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0);130#ifdef BOOST_INSTRUMENT131 std::cout << "Estimating x with Temme method 1: " << x << std::endl;132#endif133 return x;134}135//136// See:137// "Asymptotic Inversion of the Incomplete Beta Function"138// N.M. Temme139// Journal of Computation and Applied Mathematics 41 (1992) 145-157.140// Section 3.141//142template <class T, class Policy>143BOOST_MATH_GPU_ENABLED T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r, T theta, const Policy& pol)144{145 BOOST_MATH_STD_USING // ADL of std names146 147 //148 // Get first estimate for eta, see Eq 3.9 and 3.10,149 // but note there is a typo in Eq 3.10:150 //151 T eta0 = boost::math::erfc_inv(2 * z, pol);152 eta0 /= -sqrt(r / 2);153 154 T s = sin(theta);155 T c = cos(theta);156 //157 // Now we need to perturb eta0 to get eta, which we do by158 // evaluating the polynomial in 1/r at the bottom of page 151,159 // to do this we first need the error terms e1, e2 e3160 // which we'll fill into the array "terms". Since these161 // terms are themselves polynomials, we'll need another162 // array "workspace" to calculate those...163 //164 T terms[4] = { eta0 };165 T workspace[6];166 //167 // some powers of sin(theta)cos(theta) that we'll need later:168 //169 T sc = s * c;170 T sc_2 = sc * sc;171 T sc_3 = sc_2 * sc;172 T sc_4 = sc_2 * sc_2;173 T sc_5 = sc_2 * sc_3;174 T sc_6 = sc_3 * sc_3;175 T sc_7 = sc_4 * sc_3;176 //177 // Calculate e1 and put it in terms[1], see the middle of page 151:178 //179 workspace[0] = (2 * s * s - 1) / (3 * s * c);180 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co1[] = { -1, -5, 5 };181 workspace[1] = -tools::evaluate_even_polynomial(co1, s, 3) / (36 * sc_2);182 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co2[] = { 1, 21, -69, 46 };183 workspace[2] = tools::evaluate_even_polynomial(co2, s, 4) / (1620 * sc_3);184 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co3[] = { 7, -2, 33, -62, 31 };185 workspace[3] = -tools::evaluate_even_polynomial(co3, s, 5) / (6480 * sc_4);186 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co4[] = { 25, -52, -17, 88, -115, 46 };187 workspace[4] = tools::evaluate_even_polynomial(co4, s, 6) / (90720 * sc_5);188 terms[1] = tools::evaluate_polynomial(workspace, eta0, 5);189 //190 // Now evaluate e2 and put it in terms[2]:191 //192 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co5[] = { 7, 12, -78, 52 };193 workspace[0] = -tools::evaluate_even_polynomial(co5, s, 4) / (405 * sc_3);194 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co6[] = { -7, 2, 183, -370, 185 };195 workspace[1] = tools::evaluate_even_polynomial(co6, s, 5) / (2592 * sc_4);196 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co7[] = { -533, 776, -1835, 10240, -13525, 5410 };197 workspace[2] = -tools::evaluate_even_polynomial(co7, s, 6) / (204120 * sc_5);198 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co8[] = { -1579, 3747, -3372, -15821, 45588, -45213, 15071 };199 workspace[3] = -tools::evaluate_even_polynomial(co8, s, 7) / (2099520 * sc_6);200 terms[2] = tools::evaluate_polynomial(workspace, eta0, 4);201 //202 // And e3, and put it in terms[3]:203 //204 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co9[] = {449, -1259, -769, 6686, -9260, 3704 };205 workspace[0] = tools::evaluate_even_polynomial(co9, s, 6) / (102060 * sc_5);206 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co10[] = { 63149, -151557, 140052, -727469, 2239932, -2251437, 750479 };207 workspace[1] = -tools::evaluate_even_polynomial(co10, s, 7) / (20995200 * sc_6);208 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co11[] = { 29233, -78755, 105222, 146879, -1602610, 3195183, -2554139, 729754 };209 workspace[2] = tools::evaluate_even_polynomial(co11, s, 8) / (36741600 * sc_7);210 terms[3] = tools::evaluate_polynomial(workspace, eta0, 3);211 //212 // Bring the correction terms together to evaluate eta,213 // this is the last equation on page 151:214 //215 T eta = tools::evaluate_polynomial(terms, T(1/r), 4);216 //217 // Now that we have eta we need to back solve for x,218 // we seek the value of x that gives eta in Eq 3.2.219 // The two methods used are described in section 5.220 //221 // Begin by defining a few variables we'll need later:222 //223 T x;224 T s_2 = s * s;225 T c_2 = c * c;226 T alpha = c / s;227 alpha *= alpha;228 T lu = (-(eta * eta) / (2 * s_2) + log(s_2) + c_2 * log(c_2) / s_2);229 //230 // Temme doesn't specify what value to switch on here,231 // but this seems to work pretty well:232 //233 if(fabs(eta) < 0.7)234 {235 //236 // Small eta use the expansion Temme gives in the second equation237 // of section 5, it's a polynomial in eta:238 //239 workspace[0] = s * s;240 workspace[1] = s * c;241 workspace[2] = (1 - 2 * workspace[0]) / 3;242 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co12[] = { 1, -13, 13 };243 workspace[3] = tools::evaluate_polynomial(co12, workspace[0], 3) / (36 * s * c);244 static const BOOST_MATH_INT_TABLE_TYPE(T, int) co13[] = { 1, 21, -69, 46 };245 workspace[4] = tools::evaluate_polynomial(co13, workspace[0], 4) / (270 * workspace[0] * c * c);246 x = tools::evaluate_polynomial(workspace, eta, 5);247#ifdef BOOST_INSTRUMENT248 std::cout << "Estimating x with Temme method 2 (small eta): " << x << std::endl;249#endif250 }251 else252 {253 //254 // If eta is large we need to solve Eq 3.2 more directly,255 // begin by getting an initial approximation for x from256 // the last equation on page 155, this is a polynomial in u:257 //258 T u = exp(lu);259 workspace[0] = u;260 workspace[1] = alpha;261 workspace[2] = 0;262 workspace[3] = 3 * alpha * (3 * alpha + 1) / 6;263 workspace[4] = 4 * alpha * (4 * alpha + 1) * (4 * alpha + 2) / 24;264 workspace[5] = 5 * alpha * (5 * alpha + 1) * (5 * alpha + 2) * (5 * alpha + 3) / 120;265 x = tools::evaluate_polynomial(workspace, u, 6);266 //267 // At this point we may or may not have the right answer, Eq-3.2 has268 // two solutions for x for any given eta, however the mapping in 3.2269 // is 1:1 with the sign of eta and x-sin^2(theta) being the same.270 // So we can check if we have the right root of 3.2, and if not271 // switch x for 1-x. This transformation is motivated by the fact272 // that the distribution is *almost* symmetric so 1-x will be in the right273 // ball park for the solution:274 //275 if((x - s_2) * eta < 0)276 x = 1 - x;277#ifdef BOOST_INSTRUMENT278 std::cout << "Estimating x with Temme method 2 (large eta): " << x << std::endl;279#endif280 }281 //282 // The final step is a few Newton-Raphson iterations to283 // clean up our approximation for x, this is pretty cheap284 // in general, and very cheap compared to an incomplete beta285 // evaluation. The limits set on x come from the observation286 // that the sign of eta and x-sin^2(theta) are the same.287 //288 T lower, upper;289 if(eta < 0)290 {291 lower = 0;292 upper = s_2;293 }294 else295 {296 lower = s_2;297 upper = 1;298 }299 //300 // If our initial approximation is out of bounds then bisect:301 //302 if((x < lower) || (x > upper))303 x = (lower+upper) / 2;304 //305 // And iterate:306 //307#ifndef BOOST_MATH_NO_EXCEPTIONS308 try {309#endif310 x = tools::newton_raphson_iterate(311 temme_root_finder<T>(-lu, alpha), x, lower, upper, policies::digits<T, Policy>() / 2);312#ifndef BOOST_MATH_NO_EXCEPTIONS313 }314 catch (const boost::math::evaluation_error&)315 {316 // Due to numerical instability we may have cases where no root is found when317 // in fact we should just touch the origin. We simply ignore the error here318 // and return our best guess for x so far...319 // Maybe we should special case the symmetrical parameter case, but it's not clear 320 // whether that is the only situation when problems can occur.321 // See https://github.com/boostorg/math/issues/1169322 }323#endif324 return x;325}326//327// See:328// "Asymptotic Inversion of the Incomplete Beta Function"329// N.M. Temme330// Journal of Computation and Applied Mathematics 41 (1992) 145-157.331// Section 4.332//333template <class T, class Policy>334BOOST_MATH_GPU_ENABLED T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)335{336 BOOST_MATH_STD_USING // ADL of std names337 338 339 //340 // Begin by getting an initial approximation for the quantity341 // eta from the dominant part of the incomplete beta:342 //343 T eta0;344 if(p < q)345 eta0 = boost::math::gamma_q_inv(b, p, pol);346 else347 eta0 = boost::math::gamma_p_inv(b, q, pol);348 eta0 /= a;349 //350 // Define the variables and powers we'll need later on:351 //352 T mu = b / a;353 T w = sqrt(1 + mu);354 T w_2 = w * w;355 T w_3 = w_2 * w;356 T w_4 = w_2 * w_2;357 T w_5 = w_3 * w_2;358 T w_6 = w_3 * w_3;359 T w_7 = w_4 * w_3;360 T w_8 = w_4 * w_4;361 T w_9 = w_5 * w_4;362 T w_10 = w_5 * w_5;363 T d = eta0 - mu;364 T d_2 = d * d;365 T d_3 = d_2 * d;366 T d_4 = d_2 * d_2;367 T w1 = w + 1;368 T w1_2 = w1 * w1;369 T w1_3 = w1 * w1_2;370 T w1_4 = w1_2 * w1_2;371 //372 // Now we need to compute the perturbation error terms that373 // convert eta0 to eta, these are all polynomials of polynomials.374 // Probably these should be re-written to use tabulated data375 // (see examples above), but it's less of a win in this case as we376 // need to calculate the individual powers for the denominator terms377 // anyway, so we might as well use them for the numerator-polynomials378 // as well....379 //380 // Refer to p154-p155 for the details of these expansions:381 //382 T e1 = (w + 2) * (w - 1) / (3 * w);383 e1 += (w_3 + 9 * w_2 + 21 * w + 5) * d / (36 * w_2 * w1);384 e1 -= (w_4 - 13 * w_3 + 69 * w_2 + 167 * w + 46) * d_2 / (1620 * w1_2 * w_3);385 e1 -= (7 * w_5 + 21 * w_4 + 70 * w_3 + 26 * w_2 - 93 * w - 31) * d_3 / (6480 * w1_3 * w_4);386 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);387 388 T e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1) / (1620 * w1 * w_3);389 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);390 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);391 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);392 393 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);394 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);395 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);396 //397 // Combine eta0 and the error terms to compute eta (Second equation p155):398 //399 T eta = eta0 + e1 / a + e2 / (a * a) + e3 / (a * a * a);400 //401 // Now we need to solve Eq 4.2 to obtain x. For any given value of402 // eta there are two solutions to this equation, and since the distribution403 // may be very skewed, these are not related by x ~ 1-x we used when404 // implementing section 3 above. However we know that:405 //406 // cross < x <= 1 ; iff eta < mu407 // x == cross ; iff eta == mu408 // 0 <= x < cross ; iff eta > mu409 //410 // Where cross == 1 / (1 + mu)411 // Many thanks to Prof Temme for clarifying this point.412 //413 // Therefore we'll just jump straight into Newton iterations414 // to solve Eq 4.2 using these bounds, and simple bisection415 // as the first guess, in practice this converges pretty quickly416 // and we only need a few digits correct anyway:417 //418 if(eta <= 0)419 eta = tools::min_value<T>();420 T u = eta - mu * log(eta) + (1 + mu) * log(1 + mu) - mu;421 T cross = 1 / (1 + mu);422 T lower = eta < mu ? cross : 0;423 T upper = eta < mu ? 1 : cross;424 T x = (lower + upper) / 2;425 426 // Early exit for cases with numerical precision issues.427 if (cross == 0 || cross == 1) { return cross; }428 429 x = tools::newton_raphson_iterate(430 temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2);431#ifdef BOOST_INSTRUMENT432 std::cout << "Estimating x with Temme method 3: " << x << std::endl;433#endif434 return x;435}436 437template <class T, class Policy>438struct ibeta_roots439{440 BOOST_MATH_GPU_ENABLED ibeta_roots(T _a, T _b, T t, bool inv = false)441 : a(_a), b(_b), target(t), invert(inv) {}442 443 BOOST_MATH_GPU_ENABLED boost::math::tuple<T, T, T> operator()(T x)444 {445 BOOST_MATH_STD_USING // ADL of std names446 447 BOOST_FPU_EXCEPTION_GUARD448 449 T f1;450 T y = 1 - x;451 T f = ibeta_imp(a, b, x, Policy(), invert, true, &f1) - target;452 if(invert)453 f1 = -f1;454 if(y == 0)455 y = tools::min_value<T>() * 64;456 if(x == 0)457 x = tools::min_value<T>() * 64;458 459 T f2 = f1 * (-y * a + (b - 2) * x + 1);460 if(fabs(f2) < y * x * tools::max_value<T>())461 f2 /= (y * x);462 if(invert)463 f2 = -f2;464 465 // make sure we don't have a zero derivative:466 if(f1 == 0)467 f1 = (invert ? -1 : 1) * tools::min_value<T>() * 64;468 469 return boost::math::make_tuple(f, f1, f2);470 }471private:472 T a, b, target;473 bool invert;474};475 476template <class T, class Policy>477BOOST_MATH_GPU_ENABLED T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)478{479 BOOST_MATH_STD_USING // For ADL of math functions.480 481 //482 // The flag invert is set to true if we swap a for b and p for q,483 // in which case the result has to be subtracted from 1:484 //485 bool invert = false;486 //487 // Handle trivial cases first:488 //489 if(q == 0)490 {491 if(py) *py = 0;492 return 1;493 }494 else if(p == 0)495 {496 if(py) *py = 1;497 return 0;498 }499 else if(a == 1)500 {501 if(b == 1)502 {503 if(py) *py = 1 - p;504 return p;505 }506 // Change things around so we can handle as b == 1 special case below:507 BOOST_MATH_GPU_SAFE_SWAP(a, b);508 BOOST_MATH_GPU_SAFE_SWAP(p, q);509 invert = true;510 }511 //512 // Depending upon which approximation method we use, we may end up513 // calculating either x or y initially (where y = 1-x):514 //515 T x = 0; // Set to a safe zero to avoid a516 // MSVC 2005 warning C4701: potentially uninitialized local variable 'x' used517 // But code inspection appears to ensure that x IS assigned whatever the code path.518 T y;519 520 // For some of the methods we can put tighter bounds521 // on the result than simply [0,1]:522 //523 T lower = 0;524 T upper = 1;525 //526 // Student's T with b = 0.5 gets handled as a special case, swap527 // around if the arguments are in the "wrong" order:528 //529 if(a == 0.5f)530 {531 if(b == 0.5f)532 {533 x = sin(p * constants::half_pi<T>());534 x *= x;535 if(py)536 {537 *py = sin(q * constants::half_pi<T>());538 *py *= *py;539 }540 return x;541 }542 else if(b > 0.5f)543 {544 BOOST_MATH_GPU_SAFE_SWAP(a, b);545 BOOST_MATH_GPU_SAFE_SWAP(p, q);546 invert = !invert;547 }548 }549 //550 // Select calculation method for the initial estimate:551 //552 if((b == 0.5f) && (a >= 0.5f) && (p != 1))553 {554 //555 // We have a Student's T distribution:556 x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol);557 }558 else if(b == 1)559 {560 if(p < q)561 {562 if(a > 1)563 {564 x = pow(p, 1 / a);565 y = -boost::math::expm1(log(p) / a, pol);566 }567 else568 {569 x = pow(p, 1 / a);570 y = 1 - x;571 }572 }573 else574 {575 x = exp(boost::math::log1p(-q, pol) / a);576 y = -boost::math::expm1(boost::math::log1p(-q, pol) / a, pol);577 }578 if(invert)579 BOOST_MATH_GPU_SAFE_SWAP(x, y);580 if(py)581 *py = y;582 return x;583 }584 else if(a + b > 5)585 {586 //587 // When a+b is large then we can use one of Prof Temme's588 // asymptotic expansions, begin by swapping things around589 // so that p < 0.5, we do this to avoid cancellations errors590 // when p is large.591 //592 if(p > 0.5)593 {594 BOOST_MATH_GPU_SAFE_SWAP(a, b);595 BOOST_MATH_GPU_SAFE_SWAP(p, q);596 invert = !invert;597 }598 T minv = BOOST_MATH_GPU_SAFE_MIN(a, b);599 T maxv = BOOST_MATH_GPU_SAFE_MAX(a, b);600 if((sqrt(minv) > (maxv - minv)) && (minv > 5))601 {602 //603 // When a and b differ by a small amount604 // the curve is quite symmetrical and we can use an error605 // function to approximate the inverse. This is the cheapest606 // of the three Temme expansions, and the calculated value607 // for x will never be much larger than p, so we don't have608 // to worry about cancellation as long as p is small.609 //610 x = temme_method_1_ibeta_inverse(a, b, p, pol);611 y = 1 - x;612 }613 else614 {615 T r = a + b;616 T theta = asin(sqrt(a / r));617 T lambda = minv / r;618 if((lambda >= 0.2) && (lambda <= 0.8) && (r >= 10))619 {620 //621 // The second error function case is the next cheapest622 // to use, it brakes down when the result is likely to be623 // very small, if a+b is also small, but we can use a624 // cheaper expansion there in any case. As before x won't625 // be much larger than p, so as long as p is small we should626 // be free of cancellation error.627 //628 T ppa = pow(p, 1/a);629 if((ppa < 0.0025) && (a + b < 200))630 {631 x = ppa * pow(a * boost::math::beta(a, b, pol), 1/a);632 }633 else634 x = temme_method_2_ibeta_inverse(a, b, p, r, theta, pol);635 y = 1 - x;636 }637 else638 {639 //640 // If we get here then a and b are very different in magnitude641 // and we need to use the third of Temme's methods which642 // involves inverting the incomplete gamma. This is much more643 // expensive than the other methods. We also can only use this644 // method when a > b, which can lead to cancellation errors645 // if we really want y (as we will when x is close to 1), so646 // a different expansion is used in that case.647 //648 if(a < b)649 {650 BOOST_MATH_GPU_SAFE_SWAP(a, b);651 BOOST_MATH_GPU_SAFE_SWAP(p, q);652 invert = !invert;653 }654 //655 // Try and compute the easy way first:656 //657 T bet = 0;658 if (b < 2)659 {660#ifndef BOOST_MATH_NO_EXCEPTIONS661 try662#endif663 {664 bet = boost::math::beta(a, b, pol);665 666 typedef typename Policy::overflow_error_type overflow_type;667 668 BOOST_MATH_IF_CONSTEXPR(overflow_type::value != boost::math::policies::throw_on_error)669 if(bet > tools::max_value<T>())670 bet = tools::max_value<T>();671 }672#ifndef BOOST_MATH_NO_EXCEPTIONS673 catch (const std::overflow_error&)674 {675 bet = tools::max_value<T>();676 }677#endif678 }679 if(bet != 0)680 {681 y = pow(b * q * bet, 1/b);682 x = 1 - y;683 }684 else685 y = 1;686 if((y > 1e-5) && BOOST_MATH_GPU_SAFE_MIN(a, b) < 1000)687 {688 x = temme_method_3_ibeta_inverse(a, b, p, q, pol);689 y = 1 - x;690 }691 else if ((y > 1e-5) && BOOST_MATH_GPU_SAFE_MIN(a, b) > 1000)692 {693 // All options have failed, use the saddle point as a starting location:694 x = BOOST_MATH_GPU_SAFE_MAX(a, b) / (a + b);695 y = BOOST_MATH_GPU_SAFE_MIN(a, b) / (a + b);696 }697 }698 }699 }700 else if((a < 1) && (b < 1))701 {702 //703 // Both a and b less than 1,704 // there is a point of inflection at xs:705 //706 T xs = (1 - a) / (2 - a - b);707 //708 // Now we need to ensure that we start our iteration from the709 // right side of the inflection point:710 //711 T fs = boost::math::ibeta(a, b, xs, pol) - p;712 if(fabs(fs) / p < tools::epsilon<T>() * 3)713 {714 // The result is at the point of inflection, best just return it:715 *py = invert ? xs : 1 - xs;716 return invert ? 1-xs : xs;717 }718 if(fs < 0)719 {720 BOOST_MATH_GPU_SAFE_SWAP(a, b);721 BOOST_MATH_GPU_SAFE_SWAP(p, q);722 invert = !invert;723 xs = 1 - xs;724 }725 if ((a < tools::min_value<T>()) && (b > tools::min_value<T>()))726 {727 if (py)728 {729 *py = invert ? 0 : 1;730 }731 return invert ? 1 : 0; // nothing interesting going on here.732 }733 //734 // The call to beta may overflow, plus the alternative using lgamma may do the same735 // if T is a type where 1/T is infinite for small values (denorms for example).736 //737 T bet = 0;738 T xg;739 bool overflow = false;740#ifndef BOOST_MATH_NO_EXCEPTIONS741 try {742#endif743 bet = boost::math::beta(a, b, pol);744#ifndef BOOST_MATH_NO_EXCEPTIONS745 }746 catch (const std::runtime_error&)747 {748 overflow = true;749 }750#endif751 if (overflow || !(boost::math::isfinite)(bet))752 {753 xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a);754 if (xg > 2 / tools::epsilon<T>())755 xg = 2 / tools::epsilon<T>();756 }757 else758 xg = pow(a * p * bet, 1/a);759 x = xg / (1 + xg);760 y = 1 / (1 + xg);761 //762 // And finally we know that our result is below the inflection763 // point, so set an upper limit on our search:764 //765 if(x > xs)766 x = xs;767 upper = xs;768 }769 else if((a > 1) && (b > 1))770 {771 //772 // Small a and b, both greater than 1,773 // there is a point of inflection at xs,774 // and it's complement is xs2, we must always775 // start our iteration from the right side of the776 // point of inflection.777 //778 T xs = (a - 1) / (a + b - 2);779 T xs2 = (b - 1) / (a + b - 2);780 T ps = boost::math::ibeta(a, b, xs, pol) - p;781 782 if(ps < 0)783 {784 BOOST_MATH_GPU_SAFE_SWAP(a, b);785 BOOST_MATH_GPU_SAFE_SWAP(p, q);786 BOOST_MATH_GPU_SAFE_SWAP(xs, xs2);787 invert = !invert;788 }789 //790 // Estimate x and y, using expm1 to get a good estimate791 // for y when it's very small:792 //793 T lx = log(p * a * boost::math::beta(a, b, pol)) / a;794 x = exp(lx);795 y = x < 0.9 ? T(1 - x) : (T)(-boost::math::expm1(lx, pol));796 797 if((b < a) && (x < 0.2))798 {799 //800 // Under a limited range of circumstances we can improve801 // our estimate for x, frankly it's clear if this has much effect!802 //803 T ap1 = a - 1;804 T bm1 = b - 1;805 T a_2 = a * a;806 T a_3 = a * a_2;807 T b_2 = b * b;808 T terms[5] = { 0, 1 };809 terms[2] = bm1 / ap1;810 ap1 *= ap1;811 terms[3] = bm1 * (3 * a * b + 5 * b + a_2 - a - 4) / (2 * (a + 2) * ap1);812 ap1 *= (a + 1);813 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)814 / (3 * (a + 3) * (a + 2) * ap1);815 x = tools::evaluate_polynomial(terms, x, 5);816 }817 //818 // And finally we know that our result is below the inflection819 // point, so set an upper limit on our search:820 //821 if(x > xs)822 x = xs;823 upper = xs;824 }825 else /*if((a <= 1) != (b <= 1))*/826 {827 //828 // If all else fails we get here, only one of a and b829 // is above 1, and a+b is small. Start by swapping830 // things around so that we have a concave curve with b > a831 // and no points of inflection in [0,1]. As long as we expect832 // x to be small then we can use the simple (and cheap) power833 // term to estimate x, but when we expect x to be large then834 // this greatly underestimates x and leaves us trying to835 // iterate "round the corner" which may take almost forever...836 //837 // We could use Temme's inverse gamma function case in that case,838 // this works really rather well (albeit expensively) even though839 // strictly speaking we're outside it's defined range.840 //841 // However it's expensive to compute, and an alternative approach842 // which models the curve as a distorted quarter circle is much843 // cheaper to compute, and still keeps the number of iterations844 // required down to a reasonable level. With thanks to Prof Temme845 // for this suggestion.846 //847 if(b < a)848 {849 BOOST_MATH_GPU_SAFE_SWAP(a, b);850 BOOST_MATH_GPU_SAFE_SWAP(p, q);851 invert = !invert;852 }853 if (a < tools::min_value<T>())854 {855 // Avoid spurious overflows for denorms:856 if (p < 1)857 {858 x = 1;859 y = 0;860 }861 else862 {863 x = 0;864 y = 1;865 }866 }867 else if(pow(p, 1/a) < 0.5)868 {869#ifndef BOOST_MATH_NO_EXCEPTIONS870 try 871 {872#endif873 x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);874 if ((x > 1) || !(boost::math::isfinite)(x))875 x = 1;876#ifndef BOOST_MATH_NO_EXCEPTIONS877 }878 catch (const std::overflow_error&)879 {880 x = 1;881 }882#endif883 if(x == 0)884 x = boost::math::tools::min_value<T>();885 y = 1 - x;886 }887 else /*if(pow(q, 1/b) < 0.1)*/888 {889 // model a distorted quarter circle:890#ifndef BOOST_MATH_NO_EXCEPTIONS891 try 892 {893#endif894 y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);895 if ((y > 1) || !(boost::math::isfinite)(y))896 y = 1;897#ifndef BOOST_MATH_NO_EXCEPTIONS898 }899 catch (const std::overflow_error&)900 {901 y = 1;902 }903#endif904 if(y == 0)905 y = boost::math::tools::min_value<T>();906 x = 1 - y;907 }908 }909 910 //911 // Now we have a guess for x (and for y) we can set things up for912 // iteration. If x > 0.5 it pays to swap things round:913 //914 if(x > 0.5)915 {916 BOOST_MATH_GPU_SAFE_SWAP(a, b);917 BOOST_MATH_GPU_SAFE_SWAP(p, q);918 BOOST_MATH_GPU_SAFE_SWAP(x, y);919 invert = !invert;920 T l = 1 - upper;921 T u = 1 - lower;922 lower = l;923 upper = u;924 }925 //926 // lower bound for our search:927 //928 // We're not interested in denormalised answers as these tend to929 // these tend to take up lots of iterations, given that we can't get930 // accurate derivatives in this area (they tend to be infinite).931 //932 if(lower == 0)933 {934 if(invert && (py == 0))935 {936 //937 // We're not interested in answers smaller than machine epsilon:938 //939 lower = boost::math::tools::epsilon<T>();940 if(x < lower)941 x = lower;942 }943 else944 lower = boost::math::tools::min_value<T>();945 if(x < lower)946 x = lower;947 }948 boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();949 boost::math::uintmax_t max_iter_used = 0;950 //951 // Figure out how many digits to iterate towards:952 //953 int digits = boost::math::policies::digits<T, Policy>() / 2;954 if((x < 1e-50) && ((a < 1) || (b < 1)))955 {956 //957 // If we're in a region where the first derivative is very958 // large, then we have to take care that the root-finder959 // doesn't terminate prematurely. We'll bump the precision960 // up to avoid this, but we have to take care not to set the961 // precision too high or the last few iterations will just962 // thrash around and convergence may be slow in this case.963 // Try 3/4 of machine epsilon:964 //965 digits *= 3;966 digits /= 2;967 }968 //969 // Now iterate, we can use either p or q as the target here970 // depending on which is smaller:971 //972 // Since we can't use halley_iterate on device we use newton raphson973 //974 #ifndef BOOST_MATH_HAS_GPU_SUPPORT975 x = boost::math::tools::halley_iterate(976 #else977 x = boost::math::tools::newton_raphson_iterate(978 #endif979 boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);980 policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol);981 //982 // We don't really want these asserts here, but they are useful for sanity983 // checking that we have the limits right, uncomment if you suspect bugs *only*.984 //985 //BOOST_MATH_ASSERT(x != upper);986 //BOOST_MATH_ASSERT((x != lower) || (x == boost::math::tools::min_value<T>()) || (x == boost::math::tools::epsilon<T>()));987 //988 // Tidy up, if we "lower" was too high then zero is the best answer we have:989 //990 if(x == lower)991 x = 0;992 if(py)993 *py = invert ? x : 1 - x;994 return invert ? 1-x : x;995}996 997} // namespace detail998 999template <class T1, class T2, class T3, class T4, class Policy>1000BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3, T4>::type1001 ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy& pol)1002{1003 constexpr auto function = "boost::math::ibeta_inv<%1%>(%1%,%1%,%1%)";1004 BOOST_FPU_EXCEPTION_GUARD1005 typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type;1006 typedef typename policies::evaluation<result_type, Policy>::type value_type;1007 typedef typename policies::normalise<1008 Policy,1009 policies::promote_float<false>,1010 policies::promote_double<false>,1011 policies::discrete_quantile<>,1012 policies::assert_undefined<> >::type forwarding_policy;1013 1014 if(a <= 0)1015 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);1016 if(b <= 0)1017 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);1018 if((p < 0) || (p > 1))1019 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);1020 1021 value_type rx, ry;1022 1023 rx = detail::ibeta_inv_imp(1024 static_cast<value_type>(a),1025 static_cast<value_type>(b),1026 static_cast<value_type>(p),1027 static_cast<value_type>(1 - p),1028 forwarding_policy(), &ry);1029 1030 if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function);1031 return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function);1032}1033 1034template <class T1, class T2, class T3, class T4>1035BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3, T4>::type1036 ibeta_inv(T1 a, T2 b, T3 p, T4* py)1037{1038 return ibeta_inv(a, b, p, py, policies::policy<>());1039}1040 1041template <class T1, class T2, class T3>1042BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type1043 ibeta_inv(T1 a, T2 b, T3 p)1044{1045 typedef typename tools::promote_args<T1, T2, T3>::type result_type;1046 return ibeta_inv(a, b, p, static_cast<result_type*>(nullptr), policies::policy<>());1047}1048 1049template <class T1, class T2, class T3, class Policy>1050BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type1051 ibeta_inv(T1 a, T2 b, T3 p, const Policy& pol)1052{1053 typedef typename tools::promote_args<T1, T2, T3>::type result_type;1054 return ibeta_inv(a, b, p, static_cast<result_type*>(nullptr), pol);1055}1056 1057template <class T1, class T2, class T3, class T4, class Policy>1058BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3, T4>::type1059 ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy& pol)1060{1061 constexpr auto function = "boost::math::ibetac_inv<%1%>(%1%,%1%,%1%)";1062 BOOST_FPU_EXCEPTION_GUARD1063 typedef typename tools::promote_args<T1, T2, T3, T4>::type result_type;1064 typedef typename policies::evaluation<result_type, Policy>::type value_type;1065 typedef typename policies::normalise<1066 Policy,1067 policies::promote_float<false>,1068 policies::promote_double<false>,1069 policies::discrete_quantile<>,1070 policies::assert_undefined<> >::type forwarding_policy;1071 1072 if(a <= 0)1073 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);1074 if(b <= 0)1075 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);1076 if((q < 0) || (q > 1))1077 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);1078 1079 value_type rx, ry;1080 1081 rx = detail::ibeta_inv_imp(1082 static_cast<value_type>(a),1083 static_cast<value_type>(b),1084 static_cast<value_type>(1 - q),1085 static_cast<value_type>(q),1086 forwarding_policy(), &ry);1087 1088 if(py) *py = policies::checked_narrowing_cast<T4, forwarding_policy>(ry, function);1089 return policies::checked_narrowing_cast<result_type, forwarding_policy>(rx, function);1090}1091 1092template <class T1, class T2, class T3, class T4>1093BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3, T4>::type1094 ibetac_inv(T1 a, T2 b, T3 q, T4* py)1095{1096 return ibetac_inv(a, b, q, py, policies::policy<>());1097}1098 1099template <class RT1, class RT2, class RT3>1100BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type1101 ibetac_inv(RT1 a, RT2 b, RT3 q)1102{1103 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;1104 return ibetac_inv(a, b, q, static_cast<result_type*>(nullptr), policies::policy<>());1105}1106 1107template <class RT1, class RT2, class RT3, class Policy>1108BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type1109 ibetac_inv(RT1 a, RT2 b, RT3 q, const Policy& pol)1110{1111 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;1112 return ibetac_inv(a, b, q, static_cast<result_type*>(nullptr), pol);1113}1114 1115} // namespace math1116} // namespace boost1117 1118#endif // BOOST_MATH_SPECIAL_FUNCTIONS_IGAMMA_INVERSE_HPP1119 1120 1121 1122 1123