//  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, value_type q)

      : dist(d), target(p < q ? p : q), comp(p < q ? false : true) {}



   value_type operator()(value_type const& x)

   {

      return comp ? target - cdf(complement(dist, x)) : 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,

      const typename Dist::value_type& q,

      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, q);

   //

   // Max bounds of the distribution:

   //

   value_type min_bound, max_bound;

   std::tr1::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;

         }

         else

         {

            b = a;

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

            if(a < min_bound)

               a = min_bound;

            fa = f(a);

            --count;

            if(fa == 0)

               return a;

         }

      }

   }

   //

   // 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)(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)(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))

         {

            if(count == 0)

               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))

         {

            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)

               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;

   //

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

}

//

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

      const typename Dist::value_type& p,

      const typename Dist::value_type& q,

      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 <= pdf(dist, 0))

      return 0;

   return do_inverse_discrete_quantile(

      dist, 

      p, 

      q,

      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,

      const typename Dist::value_type& q,

      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

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

      return 0;

   //

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

   // upper or lower quantile:

   //

   if(p < 0.5f)

      return floor(do_inverse_discrete_quantile(

         dist, 

         p, 

         q,

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

         multiplier, 

         adder, 

         tools::equal_floor(),

         max_iter));

   // else:

   return ceil(do_inverse_discrete_quantile(

      dist, 

      p, 

      q,

      ceil(guess), 

      multiplier, 

      adder, 

      tools::equal_ceil(),

      max_iter));

}



template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      const typename Dist::value_type& q,

      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

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

      return 0;

   //

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

   // upper or lower quantile:

   //

   if(p < 0.5f)

      return ceil(do_inverse_discrete_quantile(

         dist, 

         p, 

         q,

         ceil(guess), 

         multiplier, 

         adder, 

         tools::equal_ceil(),

         max_iter));

   // else:

   return floor(do_inverse_discrete_quantile(

      dist, 

      p, 

      q,

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

      multiplier, 

      adder, 

      tools::equal_floor(),

      max_iter));

}



template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      const typename Dist::value_type& q,

      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

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

      return 0;

   return floor(do_inverse_discrete_quantile(

      dist, 

      p, 

      q,

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

      multiplier, 

      adder, 

      tools::equal_floor(),

      max_iter));

}



template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      const typename Dist::value_type& q,

      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

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

      return 0;

   return ceil(do_inverse_discrete_quantile(

      dist, 

      p, 

      q,

      ceil(guess), 

      multiplier, 

      adder, 

      tools::equal_ceil(),

      max_iter));

}



template <class Dist>

inline typename Dist::value_type 

   inverse_discrete_quantile(

      const Dist& dist,

      const typename Dist::value_type& p,

      const typename Dist::value_type& q,

      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

   if(p <= 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 floor(do_inverse_discrete_quantile(

      dist, 

      p, 

      q,

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

      multiplier, 

      adder, 

      tools::equal_nearest_integer(),

      max_iter) + 0.5f);

}



}}} // namespaces



#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_INV_DISCRETE_QUANTILE





