ceph/src/boost/libs/math/include/boost/math/distributions/detail/hypergeometric_pdf.hpp

   1 // Copyright 2008 Gautam Sewani
   2 // Copyright 2008 John Maddock
   3 //
   4 // Use, modification and distribution are subject to the
   5 // Boost Software License, Version 1.0.
   6 // (See accompanying file LICENSE_1_0.txt
   7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
  10 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
  11
  12 #include <boost/math/constants/constants.hpp>
  13 #include <boost/math/special_functions/lanczos.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15 #include <boost/math/special_functions/pow.hpp>
  16 #include <boost/math/special_functions/prime.hpp>
  17 #include <boost/math/policies/error_handling.hpp>
  18
  19 #ifdef BOOST_MATH_INSTRUMENT
  20 #include <typeinfo>
  21 #endif
  22
  23 namespace boost{ namespace math{ namespace detail{
  24
  25 template <class T, class Func>
  26 void bubble_down_one(T* first, T* last, Func f)
  27 {
  28    using std::swap;
  29    T* next = first;
  30    ++next;
  31    while((next != last) && (!f(*first, *next)))
  32    {
  33       swap(*first, *next);
  34       ++first;
  35       ++next;
  36    }
  37 }
  38
  39 template <class T>
  40 struct sort_functor
  41 {
  42    sort_functor(const T* exponents) : m_exponents(exponents){}
  43    bool operator()(int i, int j)
  44    {
  45       return m_exponents[i] > m_exponents[j];
  46    }
  47 private:
  48    const T* m_exponents;
  49 };
  50
  51 template <class T, class Lanczos, class Policy>
  52 T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const Lanczos&, const Policy&)
  53 {
  54    BOOST_MATH_STD_USING
  55
  56    BOOST_MATH_INSTRUMENT_FPU
  57    BOOST_MATH_INSTRUMENT_VARIABLE(x);
  58    BOOST_MATH_INSTRUMENT_VARIABLE(r);
  59    BOOST_MATH_INSTRUMENT_VARIABLE(n);
  60    BOOST_MATH_INSTRUMENT_VARIABLE(N);
  61    BOOST_MATH_INSTRUMENT_VARIABLE(typeid(Lanczos).name());
  62
  63    T bases[9] = {
  64       T(n) + static_cast<T>(Lanczos::g()) + 0.5f,
  65       T(r) + static_cast<T>(Lanczos::g()) + 0.5f,
  66       T(N - n) + static_cast<T>(Lanczos::g()) + 0.5f,
  67       T(N - r) + static_cast<T>(Lanczos::g()) + 0.5f,
  68       1 / (T(N) + static_cast<T>(Lanczos::g()) + 0.5f),
  69       1 / (T(x) + static_cast<T>(Lanczos::g()) + 0.5f),
  70       1 / (T(n - x) + static_cast<T>(Lanczos::g()) + 0.5f),
  71       1 / (T(r - x) + static_cast<T>(Lanczos::g()) + 0.5f),
  72       1 / (T(N - n - r + x) + static_cast<T>(Lanczos::g()) + 0.5f)
  73    };
  74    T exponents[9] = {
  75       n + T(0.5f),
  76       r + T(0.5f),
  77       N - n + T(0.5f),
  78       N - r + T(0.5f),
  79       N + T(0.5f),
  80       x + T(0.5f),
  81       n - x + T(0.5f),
  82       r - x + T(0.5f),
  83       N - n - r + x + T(0.5f)
  84    };
  85    int base_e_factors[9] = {
  86       -1, -1, -1, -1, 1, 1, 1, 1, 1
  87    };
  88    int sorted_indexes[9] = {
  89       0, 1, 2, 3, 4, 5, 6, 7, 8
  90    };
  91 #ifdef BOOST_MATH_INSTRUMENT
  92    BOOST_MATH_INSTRUMENT_FPU
  93    for(unsigned i = 0; i < 9; ++i)
  94    {
  95       BOOST_MATH_INSTRUMENT_VARIABLE(i);
  96       BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
  97       BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
  98       BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
  99       BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
 100    }
 101 #endif
 102    std::sort(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
 103 #ifdef BOOST_MATH_INSTRUMENT
 104    BOOST_MATH_INSTRUMENT_FPU
 105    for(unsigned i = 0; i < 9; ++i)
 106    {
 107       BOOST_MATH_INSTRUMENT_VARIABLE(i);
 108       BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
 109       BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
 110       BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
 111       BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
 112    }
 113 #endif
 114
 115    do{
 116       exponents[sorted_indexes[0]] -= exponents[sorted_indexes[1]];
 117       bases[sorted_indexes[1]] *= bases[sorted_indexes[0]];
 118       if((bases[sorted_indexes[1]] < tools::min_value<T>()) && (exponents[sorted_indexes[1]] != 0))
 119       {
 120          return 0;
 121       }
 122       base_e_factors[sorted_indexes[1]] += base_e_factors[sorted_indexes[0]];
 123       bubble_down_one(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
 124
 125 #ifdef BOOST_MATH_INSTRUMENT
 126       for(unsigned i = 0; i < 9; ++i)
 127       {
 128          BOOST_MATH_INSTRUMENT_VARIABLE(i);
 129          BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
 130          BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
 131          BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
 132          BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
 133       }
 134 #endif
 135    }while(exponents[sorted_indexes[1]] > 1);
 136
 137    //
 138    // Combine equal powers:
 139    //
 140    int j = 8;
 141    while(exponents[sorted_indexes[j]] == 0) --j;
 142    while(j)
 143    {
 144       while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]]))
 145       {
 146          bases[sorted_indexes[j-1]] *= bases[sorted_indexes[j]];
 147          exponents[sorted_indexes[j]] = 0;
 148          base_e_factors[sorted_indexes[j-1]] += base_e_factors[sorted_indexes[j]];
 149          bubble_down_one(sorted_indexes + j, sorted_indexes + 9, sort_functor<T>(exponents));
 150          --j;
 151       }
 152       --j;
 153
 154 #ifdef BOOST_MATH_INSTRUMENT
 155       BOOST_MATH_INSTRUMENT_VARIABLE(j);
 156       for(unsigned i = 0; i < 9; ++i)
 157       {
 158          BOOST_MATH_INSTRUMENT_VARIABLE(i);
 159          BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
 160          BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
 161          BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
 162          BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
 163       }
 164 #endif
 165    }
 166
 167 #ifdef BOOST_MATH_INSTRUMENT
 168    BOOST_MATH_INSTRUMENT_FPU
 169    for(unsigned i = 0; i < 9; ++i)
 170    {
 171       BOOST_MATH_INSTRUMENT_VARIABLE(i);
 172       BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
 173       BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
 174       BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
 175       BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
 176    }
 177 #endif
 178
 179    T result;
 180    BOOST_MATH_INSTRUMENT_VARIABLE(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])));
 181    BOOST_MATH_INSTRUMENT_VARIABLE(exponents[sorted_indexes[0]]);
 182    {
 183       BOOST_FPU_EXCEPTION_GUARD
 184       result = pow(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]);
 185    }
 186    BOOST_MATH_INSTRUMENT_VARIABLE(result);
 187    for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
 188    {
 189       BOOST_FPU_EXCEPTION_GUARD
 190       if(result < tools::min_value<T>())
 191          return 0; // short circuit further evaluation
 192       if(exponents[sorted_indexes[i]] == 1)
 193          result *= bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]]));
 194       else if(exponents[sorted_indexes[i]] == 0.5f)
 195          result *= sqrt(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])));
 196       else
 197          result *= pow(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]);
 198
 199       BOOST_MATH_INSTRUMENT_VARIABLE(result);
 200    }
 201
 202    result *= Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n + 1))
 203       * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r + 1))
 204       * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n + 1))
 205       * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - r + 1))
 206       /
 207       ( Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N + 1))
 208          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(x + 1))
 209          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n - x + 1))
 210          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r - x + 1))
 211          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n - r + x + 1)));
 212
 213    BOOST_MATH_INSTRUMENT_VARIABLE(result);
 214    return result;
 215 }
 216
 217 template <class T, class Policy>
 218 T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
 219 {
 220    BOOST_MATH_STD_USING
 221    return exp(
 222       boost::math::lgamma(T(n + 1), pol)
 223       + boost::math::lgamma(T(r + 1), pol)
 224       + boost::math::lgamma(T(N - n + 1), pol)
 225       + boost::math::lgamma(T(N - r + 1), pol)
 226       - boost::math::lgamma(T(N + 1), pol)
 227       - boost::math::lgamma(T(x + 1), pol)
 228       - boost::math::lgamma(T(n - x + 1), pol)
 229       - boost::math::lgamma(T(r - x + 1), pol)
 230       - boost::math::lgamma(T(N - n - r + x + 1), pol));
 231 }
 232
 233 template <class T>
 234 inline T integer_power(const T& x, int ex)
 235 {
 236    if(ex < 0)
 237       return 1 / integer_power(x, -ex);
 238    switch(ex)
 239    {
 240    case 0:
 241       return 1;
 242    case 1:
 243       return x;
 244    case 2:
 245       return x * x;
 246    case 3:
 247       return x * x * x;
 248    case 4:
 249       return boost::math::pow<4>(x);
 250    case 5:
 251       return boost::math::pow<5>(x);
 252    case 6:
 253       return boost::math::pow<6>(x);
 254    case 7:
 255       return boost::math::pow<7>(x);
 256    case 8:
 257       return boost::math::pow<8>(x);
 258    }
 259    BOOST_MATH_STD_USING
 260 #ifdef __SUNPRO_CC
 261    return pow(x, T(ex));
 262 #else
 263    return pow(x, ex);
 264 #endif
 265 }
 266 template <class T>
 267 struct hypergeometric_pdf_prime_loop_result_entry
 268 {
 269    T value;
 270    const hypergeometric_pdf_prime_loop_result_entry* next;
 271 };
 272
 273 #ifdef BOOST_MSVC
 274 #pragma warning(push)
 275 #pragma warning(disable:4510 4512 4610)
 276 #endif
 277
 278 struct hypergeometric_pdf_prime_loop_data
 279 {
 280    const unsigned x;
 281    const unsigned r;
 282    const unsigned n;
 283    const unsigned N;
 284    unsigned prime_index;
 285    unsigned current_prime;
 286 };
 287
 288 #ifdef BOOST_MSVC
 289 #pragma warning(pop)
 290 #endif
 291
 292 template <class T>
 293 T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry<T>& result)
 294 {
 295    while(data.current_prime <= data.N)
 296    {
 297       unsigned base = data.current_prime;
 298       int prime_powers = 0;
 299       while(base <= data.N)
 300       {
 301          prime_powers += data.n / base;
 302          prime_powers += data.r / base;
 303          prime_powers += (data.N - data.n) / base;
 304          prime_powers += (data.N - data.r) / base;
 305          prime_powers -= data.N / base;
 306          prime_powers -= data.x / base;
 307          prime_powers -= (data.n - data.x) / base;
 308          prime_powers -= (data.r - data.x) / base;
 309          prime_powers -= (data.N - data.n - data.r + data.x) / base;
 310          base *= data.current_prime;
 311       }
 312       if(prime_powers)
 313       {
 314          T p = integer_power<T>(static_cast<T>(data.current_prime), prime_powers);
 315          if((p > 1) && (tools::max_value<T>() / p < result.value))
 316          {
 317             //
 318             // The next calculation would overflow, use recursion
 319             // to sidestep the issue:
 320             //
 321             hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
 322             data.current_prime = prime(++data.prime_index);
 323             return hypergeometric_pdf_prime_loop_imp<T>(data, t);
 324          }
 325          if((p < 1) && (tools::min_value<T>() / p > result.value))
 326          {
 327             //
 328             // The next calculation would underflow, use recursion
 329             // to sidestep the issue:
 330             //
 331             hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
 332             data.current_prime = prime(++data.prime_index);
 333             return hypergeometric_pdf_prime_loop_imp<T>(data, t);
 334          }
 335          result.value *= p;
 336       }
 337       data.current_prime = prime(++data.prime_index);
 338    }
 339    //
 340    // When we get to here we have run out of prime factors,
 341    // the overall result is the product of all the partial
 342    // results we have accumulated on the stack so far, these
 343    // are in a linked list starting with "data.head" and ending
 344    // with "result".
 345    //
 346    // All that remains is to multiply them together, taking
 347    // care not to overflow or underflow.
 348    //
 349    // Enumerate partial results >= 1 in variable i
 350    // and partial results < 1 in variable j:
 351    //
 352    hypergeometric_pdf_prime_loop_result_entry<T> const *i, *j;
 353    i = &result;
 354    while(i && i->value < 1)
 355       i = i->next;
 356    j = &result;
 357    while(j && j->value >= 1)
 358       j = j->next;
 359
 360    T prod = 1;
 361
 362    while(i || j)
 363    {
 364       while(i && ((prod <= 1) || (j == 0)))
 365       {
 366          prod *= i->value;
 367          i = i->next;
 368          while(i && i->value < 1)
 369             i = i->next;
 370       }
 371       while(j && ((prod >= 1) || (i == 0)))
 372       {
 373          prod *= j->value;
 374          j = j->next;
 375          while(j && j->value >= 1)
 376             j = j->next;
 377       }
 378    }
 379
 380    return prod;
 381 }
 382
 383 template <class T, class Policy>
 384 inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
 385 {
 386    hypergeometric_pdf_prime_loop_result_entry<T> result = { 1, 0 };
 387    hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) };
 388    return hypergeometric_pdf_prime_loop_imp<T>(data, result);
 389 }
 390
 391 template <class T, class Policy>
 392 T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
 393 {
 394    BOOST_MATH_STD_USING
 395    BOOST_ASSERT(N <= boost::math::max_factorial<T>::value);
 396    T result = boost::math::unchecked_factorial<T>(n);
 397    T num[3] = {
 398       boost::math::unchecked_factorial<T>(r),
 399       boost::math::unchecked_factorial<T>(N - n),
 400       boost::math::unchecked_factorial<T>(N - r)
 401    };
 402    T denom[5] = {
 403       boost::math::unchecked_factorial<T>(N),
 404       boost::math::unchecked_factorial<T>(x),
 405       boost::math::unchecked_factorial<T>(n - x),
 406       boost::math::unchecked_factorial<T>(r - x),
 407       boost::math::unchecked_factorial<T>(N - n - r + x)
 408    };
 409    int i = 0;
 410    int j = 0;
 411    while((i < 3) || (j < 5))
 412    {
 413       while((j < 5) && ((result >= 1) || (i >= 3)))
 414       {
 415          result /= denom[j];
 416          ++j;
 417       }
 418       while((i < 3) && ((result <= 1) || (j >= 5)))
 419       {
 420          result *= num[i];
 421          ++i;
 422       }
 423    }
 424    return result;
 425 }
 426
 427
 428 template <class T, class Policy>
 429 inline typename tools::promote_args<T>::type
 430    hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
 431 {
 432    BOOST_FPU_EXCEPTION_GUARD
 433    typedef typename tools::promote_args<T>::type result_type;
 434    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 435    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
 436    typedef typename policies::normalise<
 437       Policy,
 438       policies::promote_float<false>,
 439       policies::promote_double<false>,
 440       policies::discrete_quantile<>,
 441       policies::assert_undefined<> >::type forwarding_policy;
 442
 443    value_type result;
 444    if(N <= boost::math::max_factorial<value_type>::value)
 445    {
 446       //
 447       // If N is small enough then we can evaluate the PDF via the factorials
 448       // directly: table lookup of the factorials gives the best performance
 449       // of the methods available:
 450       //
 451       result = detail::hypergeometric_pdf_factorial_imp<value_type>(x, r, n, N, forwarding_policy());
 452    }
 453    else if(N <= boost::math::prime(boost::math::max_prime - 1))
 454    {
 455       //
 456       // If N is no larger than the largest prime number in our lookup table
 457       // (104729) then we can use prime factorisation to evaluate the PDF,
 458       // this is slow but accurate:
 459       //
 460       result = detail::hypergeometric_pdf_prime_imp<value_type>(x, r, n, N, forwarding_policy());
 461    }
 462    else
 463    {
 464       //
 465       // Catch all case - use the lanczos approximation - where available -
 466       // to evaluate the ratio of factorials.  This is reasonably fast
 467       // (almost as quick as using logarithmic evaluation in terms of lgamma)
 468       // but only a few digits better in accuracy than using lgamma:
 469       //
 470       result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
 471    }
 472
 473    if(result > 1)
 474    {
 475       result = 1;
 476    }
 477    if(result < 0)
 478    {
 479       result = 0;
 480    }
 481
 482    return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)");
 483 }
 484
 485 }}} // namespaces
 486
 487 #endif
 488