text: //  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)(value_type(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)(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))
           {
              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) : (value_type)floor(guess)), 
           multiplier, 
           adder, 
           tools::equal_floor(),
           max_iter));
     // else:
     return ceil(do_inverse_discrete_quantile(
        dist, 
        p, 
        q,
        (value_type)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
  
  
  