Chris@16: //  Copyright John Maddock 2008.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
Chris@16: #define BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace detail{
Chris@16: 
Chris@16: template <class Dist>
Chris@16: struct generic_quantile_finder
Chris@16: {
Chris@16:    typedef typename Dist::value_type value_type;
Chris@16:    typedef typename Dist::policy_type policy_type;
Chris@16: 
Chris@16:    generic_quantile_finder(const Dist& d, value_type t, bool c)
Chris@16:       : dist(d), target(t), comp(c) {}
Chris@16: 
Chris@16:    value_type operator()(const value_type& x)
Chris@16:    {
Chris@16:       return comp ?
Chris@16:          value_type(target - cdf(complement(dist, x)))
Chris@16:          : value_type(cdf(dist, x) - target);
Chris@16:    }
Chris@16: 
Chris@16: private:
Chris@16:    Dist dist;
Chris@16:    value_type target;
Chris@16:    bool comp;
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T check_range_result(const T& x, const Policy& pol, const char* function)
Chris@16: {
Chris@16:    if((x >= 0) && (x < tools::min_value<T>()))
Chris@16:       return policies::raise_underflow_error<T>(function, 0, pol);
Chris@16:    if(x <= -tools::max_value<T>())
Chris@16:       return -policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:    if(x >= tools::max_value<T>())
Chris@16:       return policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class Dist>
Chris@16: typename Dist::value_type generic_quantile(const Dist& dist, const typename Dist::value_type& p, const typename Dist::value_type& guess, bool comp, const char* function)
Chris@16: {
Chris@16:    typedef typename Dist::value_type value_type;
Chris@16:    typedef typename Dist::policy_type policy_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       policy_type, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    //
Chris@16:    // Special cases first:
Chris@16:    //
Chris@16:    if(p == 0)
Chris@16:    {
Chris@16:       return comp
Chris@16:       ? check_range_result(range(dist).second, forwarding_policy(), function)
Chris@16:       : check_range_result(range(dist).first, forwarding_policy(), function);
Chris@16:    }
Chris@16:    if(p == 1)
Chris@16:    {
Chris@16:       return !comp
Chris@16:       ? check_range_result(range(dist).second, forwarding_policy(), function)
Chris@16:       : check_range_result(range(dist).first, forwarding_policy(), function);
Chris@16:    }
Chris@16: 
Chris@16:    generic_quantile_finder<Dist> f(dist, p, comp);
Chris@16:    tools::eps_tolerance<value_type> tol(policies::digits<value_type, forwarding_policy>() - 3);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_root_iterations<forwarding_policy>();
Chris@16:    std::pair<value_type, value_type> ir = tools::bracket_and_solve_root(
Chris@16:       f, guess, value_type(2), true, tol, max_iter, forwarding_policy());
Chris@16:    value_type result = ir.first + (ir.second - ir.first) / 2;
Chris@16:    if(max_iter >= policies::get_max_root_iterations<forwarding_policy>())
Chris@16:    {
Chris@101:       return policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:"
Chris@16:          " either there is no answer to quantile"
Chris@16:          " or the answer is infinite.  Current best guess is %1%", result, forwarding_policy());
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: }}} // namespaces
Chris@16: 
Chris@16: #endif // BOOST_MATH_DISTIBUTIONS_DETAIL_GENERIC_QUANTILE_HPP
Chris@16: