summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/detail/ibeta_inverse.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/special_functions/detail/ibeta_inverse.hpp')
-rw-r--r--boost/math/special_functions/detail/ibeta_inverse.hpp97
1 files changed, 73 insertions, 24 deletions
diff --git a/boost/math/special_functions/detail/ibeta_inverse.hpp b/boost/math/special_functions/detail/ibeta_inverse.hpp
index ccfa9197d9..a9fe8cd49c 100644
--- a/boost/math/special_functions/detail/ibeta_inverse.hpp
+++ b/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -35,12 +35,12 @@ struct temme_root_finder
if(y == 0)
{
T big = tools::max_value<T>() / 4;
- return boost::math::make_tuple(-big, -big);
+ return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big));
}
if(x == 0)
{
T big = tools::max_value<T>() / 4;
- return boost::math::make_tuple(-big, big);
+ return boost::math::make_tuple(static_cast<T>(-big), big);
}
T f = log(x) + a * log(y) + t;
T f1 = (1 / x) - (a / (y));
@@ -455,6 +455,11 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
BOOST_MATH_STD_USING // For ADL of math functions.
//
+ // The flag invert is set to true if we swap a for b and p for q,
+ // in which case the result has to be subtracted from 1:
+ //
+ bool invert = false;
+ //
// Handle trivial cases first:
//
if(q == 0)
@@ -467,17 +472,19 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if(py) *py = 1;
return 0;
}
- else if((a == 1) && (b == 1))
+ else if(a == 1)
{
- if(py) *py = 1 - p;
- return p;
+ if(b == 1)
+ {
+ if(py) *py = 1 - p;
+ return p;
+ }
+ // Change things around so we can handle as b == 1 special case below:
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = true;
}
//
- // The flag invert is set to true if we swap a for b and p for q,
- // in which case the result has to be subtracted from 1:
- //
- bool invert = false;
- //
// Depending upon which approximation method we use, we may end up
// calculating either x or y initially (where y = 1-x):
//
@@ -495,21 +502,61 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// Student's T with b = 0.5 gets handled as a special case, swap
// around if the arguments are in the "wrong" order:
//
- if((a == 0.5f) && (b >= 0.5f))
+ if(a == 0.5f)
{
- std::swap(a, b);
- std::swap(p, q);
- invert = !invert;
+ if(b == 0.5f)
+ {
+ x = sin(p * constants::half_pi<T>());
+ x *= x;
+ if(py)
+ {
+ *py = sin(q * constants::half_pi<T>());
+ *py *= *py;
+ }
+ return x;
+ }
+ else if(b > 0.5f)
+ {
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = !invert;
+ }
}
//
// Select calculation method for the initial estimate:
//
- if((b == 0.5f) && (a >= 0.5f))
+ if((b == 0.5f) && (a >= 0.5f) && (p != 1))
{
//
// We have a Student's T distribution:
x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol);
}
+ else if(b == 1)
+ {
+ if(p < q)
+ {
+ if(a > 1)
+ {
+ x = pow(p, 1 / a);
+ y = -boost::math::expm1(log(p) / a, pol);
+ }
+ else
+ {
+ x = pow(p, 1 / a);
+ y = 1 - x;
+ }
+ }
+ else
+ {
+ x = exp(boost::math::log1p(-q, pol) / a);
+ y = -boost::math::expm1(boost::math::log1p(-q, pol) / a, pol);
+ }
+ if(invert)
+ std::swap(x, y);
+ if(py)
+ *py = y;
+ return x;
+ }
else if(a + b > 5)
{
//
@@ -866,14 +913,16 @@ template <class T1, class T2, class T3>
inline typename tools::promote_args<T1, T2, T3>::type
ibeta_inv(T1 a, T2 b, T3 p)
{
- return ibeta_inv(a, b, p, static_cast<T1*>(0), policies::policy<>());
+ typedef typename tools::promote_args<T1, T2, T3>::type result_type;
+ return ibeta_inv(a, b, p, static_cast<result_type*>(0), policies::policy<>());
}
template <class T1, class T2, class T3, class Policy>
inline typename tools::promote_args<T1, T2, T3>::type
ibeta_inv(T1 a, T2 b, T3 p, const Policy& pol)
{
- return ibeta_inv(a, b, p, static_cast<T1*>(0), pol);
+ typedef typename tools::promote_args<T1, T2, T3>::type result_type;
+ return ibeta_inv(a, b, p, static_cast<result_type*>(0), pol);
}
template <class T1, class T2, class T3, class T4, class Policy>
@@ -892,11 +941,11 @@ inline typename tools::promote_args<T1, T2, T3, T4>::type
policies::assert_undefined<> >::type forwarding_policy;
if(a <= 0)
- 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);
+ 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);
if(b <= 0)
- 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);
+ 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);
if((q < 0) || (q > 1))
- 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);
+ 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);
value_type rx, ry;
@@ -922,16 +971,16 @@ template <class RT1, class RT2, class RT3>
inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_inv(RT1 a, RT2 b, RT3 q)
{
- typedef typename remove_cv<RT1>::type dummy;
- return ibetac_inv(a, b, q, static_cast<dummy*>(0), policies::policy<>());
+ typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
+ return ibetac_inv(a, b, q, static_cast<result_type*>(0), policies::policy<>());
}
template <class RT1, class RT2, class RT3, class Policy>
inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_inv(RT1 a, RT2 b, RT3 q, const Policy& pol)
{
- typedef typename remove_cv<RT1>::type dummy;
- return ibetac_inv(a, b, q, static_cast<dummy*>(0), pol);
+ typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
+ return ibetac_inv(a, b, q, static_cast<result_type*>(0), pol);
}
} // namespace math
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-19 00:57:59 +0000