ceph/src/boost/boost/math/special_functions/detail/ibeta_inv_ab.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 //
   7 // This is not a complete header file, it is included by beta.hpp
   8 // after it has defined it's definitions.  This inverts the incomplete
   9 // beta functions ibeta and ibetac on the first parameters "a"
  10 // and "b" using a generic root finding algorithm (TOMS Algorithm 748).
  11 //
  12
  13 #ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
  14 #define BOOST_MATH_SP_DETAIL_BETA_INV_AB
  15
  16 #ifdef _MSC_VER
  17 #pragma once
  18 #endif
  19
  20 #include <cstdint>
  21 #include <utility>
  22 #include <boost/math/tools/toms748_solve.hpp>
  23
  24 namespace boost{ namespace math{ namespace detail{
  25
  26 template <class T, class Policy>
  27 struct beta_inv_ab_t
  28 {
  29    beta_inv_ab_t(T b_, T z_, T p_, bool invert_, bool swap_ab_) : b(b_), z(z_), p(p_), invert(invert_), swap_ab(swap_ab_) {}
  30    T operator()(T a)
  31    {
  32       return invert ?
  33          p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy())
  34          : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
  35    }
  36 private:
  37    T b, z, p;
  38    bool invert, swap_ab;
  39 };
  40
  41 template <class T, class Policy>
  42 T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
  43 {
  44    BOOST_MATH_STD_USING
  45    // mean:
  46    T m = n * (sfc) / sf;
  47    T t = sqrt(n * (sfc));
  48    // standard deviation:
  49    T sigma = t / sf;
  50    // skewness
  51    T sk = (1 + sfc) / t;
  52    // kurtosis:
  53    T k = (6 - sf * (5+sfc)) / (n * (sfc));
  54    // Get the inverse of a std normal distribution:
  55    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
  56    // Set the sign:
  57    if(p < 0.5)
  58       x = -x;
  59    T x2 = x * x;
  60    // w is correction term due to skewness
  61    T w = x + sk * (x2 - 1) / 6;
  62    //
  63    // Add on correction due to kurtosis.
  64    //
  65    if(n >= 10)
  66       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
  67
  68    w = m + sigma * w;
  69    if(w < tools::min_value<T>())
  70       return tools::min_value<T>();
  71    return w;
  72 }
  73
  74 template <class T, class Policy>
  75 T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
  76 {
  77    BOOST_MATH_STD_USING  // for ADL of std lib math functions
  78    //
  79    // Special cases first:
  80    //
  81    BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
  82    if(p == 0)
  83    {
  84       return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
  85    }
  86    if(q == 0)
  87    {
  88       return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
  89    }
  90    //
  91    // Function object, this is the functor whose root
  92    // we have to solve:
  93    //
  94    beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
  95    //
  96    // Tolerance: full precision.
  97    //
  98    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
  99    //
 100    // Now figure out a starting guess for what a may be,
 101    // we'll start out with a value that'll put p or q
 102    // right bang in the middle of their range, the functions
 103    // are quite sensitive so we should need too many steps
 104    // to bracket the root from there:
 105    //
 106    T guess = 0;
 107    T factor = 5;
 108    //
 109    // Convert variables to parameters of a negative binomial distribution:
 110    //
 111    T n = b;
 112    T sf = swap_ab ? z : 1-z;
 113    T sfc = swap_ab ? 1-z : z;
 114    T u = swap_ab ? p : q;
 115    T v = swap_ab ? q : p;
 116    if(u <= pow(sf, n))
 117    {
 118       //
 119       // Result is less than 1, negative binomial approximation
 120       // is useless....
 121       //
 122       if((p < q) != swap_ab)
 123       {
 124          guess = (std::min)(T(b * 2), T(1));
 125       }
 126       else
 127       {
 128          guess = (std::min)(T(b / 2), T(1));
 129       }
 130    }
 131    if(n * n * n * u * sf > 0.005)
 132       guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);
 133
 134    if(guess < 10)
 135    {
 136       //
 137       // Negative binomial approximation not accurate in this area:
 138       //
 139       if((p < q) != swap_ab)
 140       {
 141          guess = (std::min)(T(b * 2), T(10));
 142       }
 143       else
 144       {
 145          guess = (std::min)(T(b / 2), T(10));
 146       }
 147    }
 148    else
 149       factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
 150    BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
 151    //
 152    // Max iterations permitted:
 153    //
 154    std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 155    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
 156    if(max_iter >= policies::get_max_root_iterations<Policy>())
 157       return policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
 158    return (r.first + r.second) / 2;
 159 }
 160
 161 } // namespace detail
 162
 163 template <class RT1, class RT2, class RT3, class Policy>
 164 typename tools::promote_args<RT1, RT2, RT3>::type
 165       ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
 166 {
 167    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 168    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 169    typedef typename policies::normalise<
 170       Policy,
 171       policies::promote_float<false>,
 172       policies::promote_double<false>,
 173       policies::discrete_quantile<>,
 174       policies::assert_undefined<> >::type forwarding_policy;
 175
 176    static const char* function = "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)";
 177    if(p == 0)
 178    {
 179       return policies::raise_overflow_error<result_type>(function, 0, Policy());
 180    }
 181    if(p == 1)
 182    {
 183       return tools::min_value<result_type>();
 184    }
 185
 186    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 187       detail::ibeta_inv_ab_imp(
 188          static_cast<value_type>(b),
 189          static_cast<value_type>(x),
 190          static_cast<value_type>(p),
 191          static_cast<value_type>(1 - static_cast<value_type>(p)),
 192          false, pol),
 193       function);
 194 }
 195
 196 template <class RT1, class RT2, class RT3, class Policy>
 197 typename tools::promote_args<RT1, RT2, RT3>::type
 198       ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
 199 {
 200    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 201    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 202    typedef typename policies::normalise<
 203       Policy,
 204       policies::promote_float<false>,
 205       policies::promote_double<false>,
 206       policies::discrete_quantile<>,
 207       policies::assert_undefined<> >::type forwarding_policy;
 208
 209    static const char* function = "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)";
 210    if(q == 1)
 211    {
 212       return policies::raise_overflow_error<result_type>(function, 0, Policy());
 213    }
 214    if(q == 0)
 215    {
 216       return tools::min_value<result_type>();
 217    }
 218
 219    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 220       detail::ibeta_inv_ab_imp(
 221          static_cast<value_type>(b),
 222          static_cast<value_type>(x),
 223          static_cast<value_type>(1 - static_cast<value_type>(q)),
 224          static_cast<value_type>(q),
 225          false, pol),
 226       function);
 227 }
 228
 229 template <class RT1, class RT2, class RT3, class Policy>
 230 typename tools::promote_args<RT1, RT2, RT3>::type
 231       ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
 232 {
 233    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 234    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 235    typedef typename policies::normalise<
 236       Policy,
 237       policies::promote_float<false>,
 238       policies::promote_double<false>,
 239       policies::discrete_quantile<>,
 240       policies::assert_undefined<> >::type forwarding_policy;
 241
 242    static const char* function = "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)";
 243    if(p == 0)
 244    {
 245       return tools::min_value<result_type>();
 246    }
 247    if(p == 1)
 248    {
 249       return policies::raise_overflow_error<result_type>(function, 0, Policy());
 250    }
 251
 252    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 253       detail::ibeta_inv_ab_imp(
 254          static_cast<value_type>(a),
 255          static_cast<value_type>(x),
 256          static_cast<value_type>(p),
 257          static_cast<value_type>(1 - static_cast<value_type>(p)),
 258          true, pol),
 259       function);
 260 }
 261
 262 template <class RT1, class RT2, class RT3, class Policy>
 263 typename tools::promote_args<RT1, RT2, RT3>::type
 264       ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
 265 {
 266    static const char* function = "boost::math::ibeta_invb<%1%>(%1%, %1%, %1%)";
 267    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 268    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 269    typedef typename policies::normalise<
 270       Policy,
 271       policies::promote_float<false>,
 272       policies::promote_double<false>,
 273       policies::discrete_quantile<>,
 274       policies::assert_undefined<> >::type forwarding_policy;
 275
 276    if(q == 1)
 277    {
 278       return tools::min_value<result_type>();
 279    }
 280    if(q == 0)
 281    {
 282       return policies::raise_overflow_error<result_type>(function, 0, Policy());
 283    }
 284
 285    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 286       detail::ibeta_inv_ab_imp(
 287          static_cast<value_type>(a),
 288          static_cast<value_type>(x),
 289          static_cast<value_type>(1 - static_cast<value_type>(q)),
 290          static_cast<value_type>(q),
 291          true, pol),
 292          function);
 293 }
 294
 295 template <class RT1, class RT2, class RT3>
 296 inline typename tools::promote_args<RT1, RT2, RT3>::type
 297          ibeta_inva(RT1 b, RT2 x, RT3 p)
 298 {
 299    return boost::math::ibeta_inva(b, x, p, policies::policy<>());
 300 }
 301
 302 template <class RT1, class RT2, class RT3>
 303 inline typename tools::promote_args<RT1, RT2, RT3>::type
 304          ibetac_inva(RT1 b, RT2 x, RT3 q)
 305 {
 306    return boost::math::ibetac_inva(b, x, q, policies::policy<>());
 307 }
 308
 309 template <class RT1, class RT2, class RT3>
 310 inline typename tools::promote_args<RT1, RT2, RT3>::type
 311          ibeta_invb(RT1 a, RT2 x, RT3 p)
 312 {
 313    return boost::math::ibeta_invb(a, x, p, policies::policy<>());
 314 }
 315
 316 template <class RT1, class RT2, class RT3>
 317 inline typename tools::promote_args<RT1, RT2, RT3>::type
 318          ibetac_invb(RT1 a, RT2 x, RT3 q)
 319 {
 320    return boost::math::ibetac_invb(a, x, q, policies::policy<>());
 321 }
 322
 323 } // namespace math
 324 } // namespace boost
 325
 326 #endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB
 327
 328
 329