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