SV platform dependency builds

//  Copyright John Maddock 2007.

//  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_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE

#define BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE

#include <algorithm>

namespace boost{ namespace math{ namespace detail{

//

// Functor for root finding algorithm:

//

template <class Dist>

struct distribution_quantile_finder

{

   typedef typename Dist::value_type value_type;

   typedef typename Dist::policy_type policy_type;

   distribution_quantile_finder(const Dist d, value_type p, bool c)

      : dist(d), target(p), comp(c) {}

   value_type operator()(value_type const& x)

   {

      return comp ? value_type(target - cdf(complement(dist, x))) : value_type(cdf(dist, x) - target);

   }

private:

   Dist dist;

   value_type target;

   bool comp;

};

//

// The purpose of adjust_bounds, is to toggle the last bit of the

// range so that both ends round to the same integer, if possible.

// If they do both round the same then we terminate the search

// for the root *very* quickly when finding an integer result.

// At the point that this function is called we know that "a" is

// below the root and "b" above it, so this change can not result

// in the root no longer being bracketed.

//

template <class Real, class Tol>

void adjust_bounds(Real& /* a */, Real& /* b */, Tol const& /* tol */){}

template <class Real>

void adjust_bounds(Real& /* a */, Real& b, tools::equal_floor const& /* tol */)

{

   BOOST_MATH_STD_USING

   b -= tools::epsilon<Real>() * b;

}

template <class Real>

void adjust_bounds(Real& a, Real& /* b */, tools::equal_ceil const& /* tol */)

{

   BOOST_MATH_STD_USING

   a += tools::epsilon<Real>() * a;

}

template <class Real>

void adjust_bounds(Real& a, Real& b, tools::equal_nearest_integer const& /* tol */)

{

   BOOST_MATH_STD_USING

   a += tools::epsilon<Real>() * a;

   b -= tools::epsilon<Real>() * b;

}

//

// This is where all the work is done:

//

template <class Dist, class Tolerance>

typename Dist::value_type 

   do_inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      bool comp,

      typename Dist::value_type guess,

      const typename Dist::value_type& multiplier,

      typename Dist::value_type adder,

      const Tolerance& tol,

      boost::uintmax_t& max_iter)

{

   typedef typename Dist::value_type value_type;

   typedef typename Dist::policy_type policy_type;

   static const char* function = "boost::math::do_inverse_discrete_quantile<%1%>";

   BOOST_MATH_STD_USING

   distribution_quantile_finder<Dist> f(dist, p, comp);

   //

   // Max bounds of the distribution:

   //

   value_type min_bound, max_bound;

   boost::math::tie(min_bound, max_bound) = support(dist);

   if(guess > max_bound)

      guess = max_bound;

   if(guess < min_bound)

      guess = min_bound;

   value_type fa = f(guess);

   boost::uintmax_t count = max_iter - 1;

   value_type fb(fa), a(guess), b =0; // Compiler warning C4701: potentially uninitialized local variable 'b' used

   if(fa == 0)

      return guess;

   //

   // For small expected results, just use a linear search:

   //

   if(guess < 10)

   {

      b = a;

      while((a < 10) && (fa * fb >= 0))

      {

         if(fb <= 0)

         {

            a = b;

            b = a + 1;

            if(b > max_bound)

               b = max_bound;

            fb = f(b);

            --count;

            if(fb == 0)

               return b;

            if(a == b)

               return b; // can't go any higher!

         }

         else

         {

            b = a;

            a = (std::max)(value_type(b - 1), value_type(0));

            if(a < min_bound)

               a = min_bound;

            fa = f(a);

            --count;

            if(fa == 0)

               return a;

            if(a == b)

               return a;  //  We can't go any lower than this!

         }

      }

   }

   //

   // Try and bracket using a couple of additions first, 

   // we're assuming that "guess" is likely to be accurate

   // to the nearest int or so:

   //

   else if(adder != 0)

   {

      //

      // If we're looking for a large result, then bump "adder" up

      // by a bit to increase our chances of bracketing the root:

      //

      //adder = (std::max)(adder, 0.001f * guess);

      if(fa < 0)

      {

         b = a + adder;

         if(b > max_bound)

            b = max_bound;

      }

      else

      {

         b = (std::max)(value_type(a - adder), value_type(0));

         if(b < min_bound)

            b = min_bound;

      }

      fb = f(b);

      --count;

      if(fb == 0)

         return b;

      if(count && (fa * fb >= 0))

      {

         //

         // We didn't bracket the root, try 

         // once more:

         //

         a = b;

         fa = fb;

         if(fa < 0)

         {

            b = a + adder;

            if(b > max_bound)

               b = max_bound;

         }

         else

         {

            b = (std::max)(value_type(a - adder), value_type(0));

            if(b < min_bound)

               b = min_bound;

         }

         fb = f(b);

         --count;

      }

      if(a > b)

      {

         using std::swap;

         swap(a, b);

         swap(fa, fb);

      }

   }

   //

   // If the root hasn't been bracketed yet, try again

   // using the multiplier this time:

   //

   if((boost::math::sign)(fb) == (boost::math::sign)(fa))

   {

      if(fa < 0)

      {

         //

         // Zero is to the right of x2, so walk upwards

         // until we find it:

         //

         while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))

         {

            if(count == 0)

               return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type());

            a = b;

            fa = fb;

            b *= multiplier;

            if(b > max_bound)

               b = max_bound;

            fb = f(b);

            --count;

            BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);

         }

      }

      else

      {

         //

         // Zero is to the left of a, so walk downwards

         // until we find it:

         //

         while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b))

         {

            if(fabs(a) < tools::min_value<value_type>())

            {

               // Escape route just in case the answer is zero!

               max_iter -= count;

               max_iter += 1;

               return 0;

            }

            if(count == 0)

               return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type());

            b = a;

            fb = fa;

            a /= multiplier;

            if(a < min_bound)

               a = min_bound;

            fa = f(a);

            --count;

            BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);

         }

      }

   }

   max_iter -= count;

   if(fa == 0)

      return a;

   if(fb == 0)

      return b;

   if(a == b)

      return b;  // Ran out of bounds trying to bracket - there is no answer!

   //

   // Adjust bounds so that if we're looking for an integer

   // result, then both ends round the same way:

   //

   adjust_bounds(a, b, tol);

   //

   // We don't want zero or denorm lower bounds:

   //

   if(a < tools::min_value<value_type>())

      a = tools::min_value<value_type>();

   //

   // Go ahead and find the root:

   //

   std::pair<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());

   max_iter += count;

   BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);

   return (r.first + r.second) / 2;

}

//

// Some special routine for rounding up and down:

// We want to check and see if we are very close to an integer, and if so test to see if

// that integer is an exact root of the cdf.  We do this because our root finder only

// guarantees to find *a root*, and there can sometimes be many consecutive floating

// point values which are all roots.  This is especially true if the target probability

// is very close 1.

//

template <class Dist>

inline typename Dist::value_type round_to_floor(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)

{

   BOOST_MATH_STD_USING

   typename Dist::value_type cc = ceil(result);

   typename Dist::value_type pp = cc <= support(d).second ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 1;

   if(pp == p)

      result = cc;

   else

      result = floor(result);

   //

   // Now find the smallest integer <= result for which we get an exact root:

   //

   while(result != 0)

   {

      cc = result - 1;

      if(cc < support(d).first)

         break;

      pp = c ? cdf(complement(d, cc)) : cdf(d, cc);

      if(pp == p)

         result = cc;

      else if(c ? pp > p : pp < p)

         break;

      result -= 1;

   }

   return result;

}

#ifdef BOOST_MSVC

#pragma warning(push)

#pragma warning(disable:4127)

#endif

template <class Dist>

inline typename Dist::value_type round_to_ceil(const Dist& d, typename Dist::value_type result, typename Dist::value_type p, bool c)

{

   BOOST_MATH_STD_USING

   typename Dist::value_type cc = floor(result);

   typename Dist::value_type pp = cc >= support(d).first ? c ? cdf(complement(d, cc)) : cdf(d, cc) : 0;

   if(pp == p)

      result = cc;

   else

      result = ceil(result);

   //

   // Now find the largest integer >= result for which we get an exact root:

   //

   while(true)

   {

      cc = result + 1;

      if(cc > support(d).second)

         break;

      pp = c ? cdf(complement(d, cc)) : cdf(d, cc);

      if(pp == p)

         result = cc;

      else if(c ? pp < p : pp > p)

         break;

      result += 1;

   }

   return result;

}

#ifdef BOOST_MSVC

#pragma warning(pop)

#endif

//

// Now finally are the public API functions.

// There is one overload for each policy,

// each one is responsible for selecting the correct

// termination condition, and rounding the result

// to an int where required.

//

template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      typename Dist::value_type p,

      bool c,

      const typename Dist::value_type& guess,

      const typename Dist::value_type& multiplier,

      const typename Dist::value_type& adder,

      const policies::discrete_quantile<policies::real>&,

      boost::uintmax_t& max_iter)

{

   if(p > 0.5)

   {

      p = 1 - p;

      c = !c;

   }

   typename Dist::value_type pp = c ? 1 - p : p;

   if(pp <= pdf(dist, 0))

      return 0;

   return do_inverse_discrete_quantile(

      dist, 

      p, 

      c,

      guess, 

      multiplier, 

      adder, 

      tools::eps_tolerance<typename Dist::value_type>(policies::digits<typename Dist::value_type, typename Dist::policy_type>()),

      max_iter);

}

template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      bool c,

      const typename Dist::value_type& guess,

      const typename Dist::value_type& multiplier,

      const typename Dist::value_type& adder,

      const policies::discrete_quantile<policies::integer_round_outwards>&,

      boost::uintmax_t& max_iter)

{

   typedef typename Dist::value_type value_type;

   BOOST_MATH_STD_USING

   typename Dist::value_type pp = c ? 1 - p : p;

   if(pp <= pdf(dist, 0))

      return 0;

   //

   // What happens next depends on whether we're looking for an 

   // upper or lower quantile:

   //

   if(pp < 0.5f)

      return round_to_floor(dist, do_inverse_discrete_quantile(

         dist, 

         p, 

         c,

         (guess < 1 ? value_type(1) : (value_type)floor(guess)), 

         multiplier, 

         adder, 

         tools::equal_floor(),

         max_iter), p, c);

   // else:

   return round_to_ceil(dist, do_inverse_discrete_quantile(

      dist, 

      p, 

      c,

      (value_type)ceil(guess), 

      multiplier, 

      adder, 

      tools::equal_ceil(),

      max_iter), p, c);

}

template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      bool c,

      const typename Dist::value_type& guess,

      const typename Dist::value_type& multiplier,

      const typename Dist::value_type& adder,

      const policies::discrete_quantile<policies::integer_round_inwards>&,

      boost::uintmax_t& max_iter)

{

   typedef typename Dist::value_type value_type;

   BOOST_MATH_STD_USING

   typename Dist::value_type pp = c ? 1 - p : p;

   if(pp <= pdf(dist, 0))

      return 0;

   //

   // What happens next depends on whether we're looking for an 

   // upper or lower quantile:

   //

   if(pp < 0.5f)

      return round_to_ceil(dist, do_inverse_discrete_quantile(

         dist, 

         p, 

         c,

         ceil(guess), 

         multiplier, 

         adder, 

         tools::equal_ceil(),

         max_iter), p, c);

   // else:

   return round_to_floor(dist, do_inverse_discrete_quantile(

      dist, 

      p, 

      c,

      (guess < 1 ? value_type(1) : floor(guess)), 

      multiplier, 

      adder, 

      tools::equal_floor(),

      max_iter), p, c);

}

template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      bool c,

      const typename Dist::value_type& guess,

      const typename Dist::value_type& multiplier,

      const typename Dist::value_type& adder,

      const policies::discrete_quantile<policies::integer_round_down>&,

      boost::uintmax_t& max_iter)

{

   typedef typename Dist::value_type value_type;

   BOOST_MATH_STD_USING

   typename Dist::value_type pp = c ? 1 - p : p;

   if(pp <= pdf(dist, 0))

      return 0;

   return round_to_floor(dist, do_inverse_discrete_quantile(

      dist, 

      p, 

      c,

      (guess < 1 ? value_type(1) : floor(guess)), 

      multiplier, 

      adder, 

      tools::equal_floor(),

      max_iter), p, c);

}

template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      bool c,

      const typename Dist::value_type& guess,

      const typename Dist::value_type& multiplier,

      const typename Dist::value_type& adder,

      const policies::discrete_quantile<policies::integer_round_up>&,

      boost::uintmax_t& max_iter)

{

   BOOST_MATH_STD_USING

   typename Dist::value_type pp = c ? 1 - p : p;

   if(pp <= pdf(dist, 0))

      return 0;

   return round_to_ceil(dist, do_inverse_discrete_quantile(

      dist, 

      p, 

      c,

      ceil(guess), 

      multiplier, 

      adder, 

      tools::equal_ceil(),

      max_iter), p, c);

}

template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      bool c,

      const typename Dist::value_type& guess,

      const typename Dist::value_type& multiplier,

      const typename Dist::value_type& adder,

      const policies::discrete_quantile<policies::integer_round_nearest>&,

      boost::uintmax_t& max_iter)

{

   typedef typename Dist::value_type value_type;

   BOOST_MATH_STD_USING

   typename Dist::value_type pp = c ? 1 - p : p;

   if(pp <= pdf(dist, 0))

      return 0;

   //

   // Note that we adjust the guess to the nearest half-integer:

   // this increase the chances that we will bracket the root

   // with two results that both round to the same integer quickly.

   //

   return round_to_floor(dist, do_inverse_discrete_quantile(

      dist, 

      p, 

      c,

      (guess < 0.5f ? value_type(1.5f) : floor(guess + 0.5f) + 0.5f), 

      multiplier, 

      adder, 

      tools::equal_nearest_integer(),

      max_iter) + 0.5f, p, c);

}

}}} // namespaces

#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE