text: // boost\math\distributions\non_central_chi_squared.hpp
  
  // Copyright John Maddock 2008.
  
  // 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_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
  #define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
  
  #include <boost/math/distributions/fwd.hpp>
  #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
  #include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
  #include <boost/math/special_functions/round.hpp> // for iround
  #include <boost/math/distributions/complement.hpp> // complements
  #include <boost/math/distributions/chi_squared.hpp> // central distribution
  #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
  #include <boost/math/special_functions/fpclassify.hpp> // isnan.
  #include <boost/math/tools/roots.hpp> // for root finding.
  #include <boost/math/distributions/detail/generic_mode.hpp>
  #include <boost/math/distributions/detail/generic_quantile.hpp>
  
  namespace boost
  {
     namespace math
     {
  
        template <class RealType, class Policy>
        class non_central_chi_squared_distribution;
  
        namespace detail{
  
           template <class T, class Policy>
           T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0)
           {
              //
              // Computes the complement of the Non-Central Chi-Square
              // Distribution CDF by summing a weighted sum of complements
              // of the central-distributions.  The weighting factor is
              // a Poisson Distribution.
              //
              // This is an application of the technique described in:
              //
              // Computing discrete mixtures of continuous
              // distributions: noncentral chisquare, noncentral t
              // and the distribution of the square of the sample
              // multiple correlation coeficient.
              // D. Benton, K. Krishnamoorthy.
              // Computational Statistics & Data Analysis 43 (2003) 249 - 267
              //
              BOOST_MATH_STD_USING
  
              // Special case:
              if(x == 0)
                 return 1;
  
              //
              // Initialize the variables we'll be using:
              //
              T lambda = theta / 2;
              T del = f / 2;
              T y = x / 2;
              boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
              T errtol = boost::math::policies::get_epsilon<T, Policy>();
              T sum = init_sum;
              //
              // k is the starting location for iteration, we'll
              // move both forwards and backwards from this point.
              // k is chosen as the peek of the Poisson weights, which
              // will occur *before* the largest term.
              //
              int k = iround(lambda, pol);
              // Forwards and backwards Poisson weights:
              T poisf = boost::math::gamma_p_derivative(1 + k, lambda, pol);
              T poisb = poisf * k / lambda;
              // Initial forwards central chi squared term:
              T gamf = boost::math::gamma_q(del + k, y, pol);
              // Forwards and backwards recursion terms on the central chi squared:
              T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
              T xtermb = xtermf * (del + k) / y;
              // Initial backwards central chi squared term:
              T gamb = gamf - xtermb;
  
              //
              // Forwards iteration first, this is the
              // stable direction for the gamma function
              // recurrences:
              //
              int i;
              for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i)
              {
                 T term = poisf * gamf;
                 sum += term;
                 poisf *= lambda / (i + 1);
                 gamf += xtermf;
                 xtermf *= y / (del + i + 1);
                 if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf))
                    break;
              }
              //Error check:
              if(static_cast<boost::uintmax_t>(i-k) >= max_iter)
                 policies::raise_evaluation_error(
                    "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
                    "Series did not converge, closest value was %1%", sum, pol);
              //
              // Now backwards iteration: the gamma
              // function recurrences are unstable in this
              // direction, we rely on the terms deminishing in size
              // faster than we introduce cancellation errors.
              // For this reason it's very important that we start
              // *before* the largest term so that backwards iteration
              // is strictly converging.
              //
              for(i = k - 1; i >= 0; --i)
              {
                 T term = poisb * gamb;
                 sum += term;
                 poisb *= i / lambda;
                 xtermb *= (del + i) / y;
                 gamb -= xtermb;
                 if((sum == 0) || (fabs(term / sum) < errtol))
                    break;
              }
  
              return sum;
           }
  
           template <class T, class Policy>
           T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0)
           {
              //
              // This is an implementation of:
              //
              // Algorithm AS 275:
              // Computing the Non-Central #2 Distribution Function
              // Cherng G. Ding
              // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
              //
              // This uses a stable forward iteration to sum the
              // CDF, unfortunately this can not be used for large
              // values of the non-centrality parameter because:
              // * The first term may underfow to zero.
              // * We may need an extra-ordinary number of terms
              //   before we reach the first *significant* term.
              //
              BOOST_MATH_STD_USING
              // Special case:
              if(x == 0)
                 return 0;
              T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
              T lambda = theta / 2;
              T vk = exp(-lambda);
              T uk = vk;
              T sum = init_sum + tk * vk;
              if(sum == 0)
                 return sum;
  
              boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
              T errtol = boost::math::policies::get_epsilon<T, Policy>();
  
              int i;
              T lterm(0), term(0);
              for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
              {
                 tk = tk * x / (f + 2 * i);
                 uk = uk * lambda / i;
                 vk = vk + uk;
                 lterm = term;
                 term = vk * tk;
                 sum += term;
                 if((fabs(term / sum) < errtol) && (term <= lterm))
                    break;
              }
              //Error check:
              if(static_cast<boost::uintmax_t>(i) >= max_iter)
                 policies::raise_evaluation_error(
                    "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
                    "Series did not converge, closest value was %1%", sum, pol);
              return sum;
           }
  
  
           template <class T, class Policy>
           T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum)
           {
              //
              // This is taken more or less directly from:
              //
              // Computing discrete mixtures of continuous
              // distributions: noncentral chisquare, noncentral t
              // and the distribution of the square of the sample
              // multiple correlation coeficient.
              // D. Benton, K. Krishnamoorthy.
              // Computational Statistics & Data Analysis 43 (2003) 249 - 267
              //
              // We're summing a Poisson weighting term multiplied by
              // a central chi squared distribution.
              //
              BOOST_MATH_STD_USING
              // Special case:
              if(y == 0)
                 return 0;
              boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
              T errtol = boost::math::policies::get_epsilon<T, Policy>();
              T errorf(0), errorb(0);
  
              T x = y / 2;
              T del = lambda / 2;
              //
              // Starting location for the iteration, we'll iterate
              // both forwards and backwards from this point.  The
              // location chosen is the maximum of the Poisson weight
              // function, which ocurrs *after* the largest term in the
              // sum.
              //
              int k = iround(del, pol);
              T a = n / 2 + k;
              // Central chi squared term for forward iteration:
              T gamkf = boost::math::gamma_p(a, x, pol);
  
              if(lambda == 0)
                 return gamkf;
              // Central chi squared term for backward iteration:
              T gamkb = gamkf;
              // Forwards Poisson weight:
              T poiskf = gamma_p_derivative(k+1, del, pol);
              // Backwards Poisson weight:
              T poiskb = poiskf;
              // Forwards gamma function recursion term:
              T xtermf = boost::math::gamma_p_derivative(a, x, pol);
              // Backwards gamma function recursion term:
              T xtermb = xtermf * x / a;
              T sum = init_sum + poiskf * gamkf;
              if(sum == 0)
                 return sum;
              int i = 1;
              //
              // Backwards recursion first, this is the stable
              // direction for gamma function recurrences:
              //
              while(i <= k)
              {
                 xtermb *= (a - i + 1) / x;
                 gamkb += xtermb;
                 poiskb = poiskb * (k - i + 1) / del;
                 errorf = errorb;
                 errorb = gamkb * poiskb;
                 sum += errorb;
                 if((fabs(errorb / sum) < errtol) && (errorb <= errorf))
                    break;
                 ++i;
              }
              i = 1;
              //
              // Now forwards recursion, the gamma function
              // recurrence relation is unstable in this direction,
              // so we rely on the magnitude of successive terms
              // decreasing faster than we introduce cancellation error.
              // For this reason it's vital that k is chosen to be *after*
              // the largest term, so that successive forward iterations
              // are strictly (and rapidly) converging.
              //
              do
              {
                 xtermf = xtermf * x / (a + i - 1);
                 gamkf = gamkf - xtermf;
                 poiskf = poiskf * del / (k + i);
                 errorf = poiskf * gamkf;
                 sum += errorf;
                 ++i;
              }while((fabs(errorf / sum) > errtol) && (static_cast<boost::uintmax_t>(i) < max_iter));
  
              //Error check:
              if(static_cast<boost::uintmax_t>(i) >= max_iter)
                 policies::raise_evaluation_error(
                    "cdf(non_central_chi_squared_distribution<%1%>, %1%)",
                    "Series did not converge, closest value was %1%", sum, pol);
  
              return sum;
           }
  
           template <class T, class Policy>
           T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
           {
              //
              // As above but for the PDF:
              //
              BOOST_MATH_STD_USING
              boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
              T errtol = boost::math::policies::get_epsilon<T, Policy>();
              T x2 = x / 2;
              T n2 = n / 2;
              T l2 = lambda / 2;
              T sum = 0;
              int k = itrunc(l2);
              T pois = gamma_p_derivative(k + 1, l2, pol) * gamma_p_derivative(n2 + k, x2);
              if(pois == 0)
                 return 0;
              T poisb = pois;
              for(int i = k; ; ++i)
              {
                 sum += pois;
                 if(pois / sum < errtol)
                    break;
                 if(static_cast<boost::uintmax_t>(i - k) >= max_iter)
                    return policies::raise_evaluation_error(
                       "pdf(non_central_chi_squared_distribution<%1%>, %1%)",
                       "Series did not converge, closest value was %1%", sum, pol);
                 pois *= l2 * x2 / ((i + 1) * (n2 + i));
              }
              for(int i = k - 1; i >= 0; --i)
              {
                 poisb *= (i + 1) * (n2 + i) / (l2 * x2);
                 sum += poisb;
                 if(poisb / sum < errtol)
                    break;
              }
              return sum / 2;
           }
  
           template <class RealType, class Policy>
           inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
           {
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
  
              BOOST_MATH_STD_USING
              value_type result;
              if(l == 0)
                 result = cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x);
              else if(x > k + l)
              {
                 // Complement is the smaller of the two:
                 result = detail::non_central_chi_square_q(
                    static_cast<value_type>(x),
                    static_cast<value_type>(k),
                    static_cast<value_type>(l),
                    forwarding_policy(),
                    static_cast<value_type>(invert ? 0 : -1));
                 invert = !invert;
              }
              else if(l < 200)
              {
                 // For small values of the non-centrality parameter
                 // we can use Ding's method:
                 result = detail::non_central_chi_square_p_ding(
                    static_cast<value_type>(x),
                    static_cast<value_type>(k),
                    static_cast<value_type>(l),
                    forwarding_policy(),
                    static_cast<value_type>(invert ? -1 : 0));
              }
              else
              {
                 // For largers values of the non-centrality
                 // parameter Ding's method will consume an
                 // extra-ordinary number of terms, and worse
                 // may return zero when the result is in fact
                 // finite, use Krishnamoorthy's method instead:
                 result = detail::non_central_chi_square_p(
                    static_cast<value_type>(x),
                    static_cast<value_type>(k),
                    static_cast<value_type>(l),
                    forwarding_policy(),
                    static_cast<value_type>(invert ? -1 : 0));
              }
              if(invert)
                 result = -result;
              return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 result,
                 "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
           }
  
           template <class T, class Policy>
           struct nccs_quantile_functor
           {
              nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
                 : dist(d), target(t), comp(c) {}
  
              T operator()(const T& x)
              {
                 return comp ?
                    target - cdf(complement(dist, x))
                    : cdf(dist, x) - target;
              }
  
           private:
              non_central_chi_squared_distribution<T,Policy> dist;
              T target;
              bool comp;
           };
  
           template <class RealType, class Policy>
           RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
           {
              static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
  
              value_type k = dist.degrees_of_freedom();
              value_type l = dist.non_centrality();
              value_type r;
              if(!detail::check_df(
                 function,
                 k, &r, Policy())
                 ||
              !detail::check_non_centrality(
                 function,
                 l,
                 &r,
                 Policy())
                 ||
              !detail::check_probability(
                 function,
                 static_cast<value_type>(p),
                 &r,
                 Policy()))
                    return (RealType)r;
  
              value_type b = (l * l) / (k + 3 * l);
              value_type c = (k + 3 * l) / (k + 2 * l);
              value_type ff = (k + 2 * l) / (c * c);
              value_type guess;
              if(comp)
                 guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
              else
                 guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
  
              if(guess < 0)
                 guess = tools::min_value<value_type>();
  
              value_type result = detail::generic_quantile(
                 non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),
                 p,
                 guess,
                 comp,
                 function);
              return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 result,
                 function);
           }
  
           template <class RealType, class Policy>
           RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
           {
              BOOST_MATH_STD_USING
              static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)";
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
  
              value_type k = dist.degrees_of_freedom();
              value_type l = dist.non_centrality();
              value_type r;
              if(!detail::check_df(
                 function,
                 k, &r, Policy())
                 ||
              !detail::check_non_centrality(
                 function,
                 l,
                 &r,
                 Policy())
                 ||
              !detail::check_positive_x(
                 function,
                 (value_type)x,
                 &r,
                 Policy()))
                    return (RealType)r;
  
           if(l == 0)
              return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
  
           // Special case:
           if(x == 0)
              return 0;
           if(l > 50)
           {
              r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
           }
           else
           {
              r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
              if(fabs(r) >= tools::log_max_value<RealType>() / 4)
              {
                 r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
              }
              else
              {
                 r = exp(r);
                 r = 0.5f * r
                    * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
              }
           }
           return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 r,
                 function);
           }
  
           template <class RealType, class Policy>
           struct degrees_of_freedom_finder
           {
              degrees_of_freedom_finder(
                 RealType lam_, RealType x_, RealType p_, bool c)
                 : lam(lam_), x(x_), p(p_), comp(c) {}
  
              RealType operator()(const RealType& v)
              {
                 non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
                 return comp ?
                    p - cdf(complement(d, x))
                    : cdf(d, x) - p;
              }
           private:
              RealType lam;
              RealType x;
              RealType p;
              bool comp;
           };
  
           template <class RealType, class Policy>
           inline RealType find_degrees_of_freedom(
              RealType lam, RealType x, RealType p, RealType q, const Policy& pol)
           {
              const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
              if((p == 0) || (q == 0))
              {
                 //
                 // Can't a thing if one of p and q is zero:
                 //
                 return policies::raise_evaluation_error<RealType>(function,
                    "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
                    RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
              }
              degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
              tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
              boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
              //
              // Pick an initial guess that we know will give us a probability
              // right around 0.5.
              //
              RealType guess = x - lam;
              if(guess < 1)
                 guess = 1;
              std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
                 f, guess, RealType(2), false, tol, max_iter, pol);
              RealType result = ir.first + (ir.second - ir.first) / 2;
              if(max_iter >= policies::get_max_root_iterations<Policy>())
              {
                 policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
                    " or there is no answer to problem.  Current best guess is %1%", result, Policy());
              }
              return result;
           }
  
           template <class RealType, class Policy>
           struct non_centrality_finder
           {
              non_centrality_finder(
                 RealType v_, RealType x_, RealType p_, bool c)
                 : v(v_), x(x_), p(p_), comp(c) {}
  
              RealType operator()(const RealType& lam)
              {
                 non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
                 return comp ?
                    p - cdf(complement(d, x))
                    : cdf(d, x) - p;
              }
           private:
              RealType v;
              RealType x;
              RealType p;
              bool comp;
           };
  
           template <class RealType, class Policy>
           inline RealType find_non_centrality(
              RealType v, RealType x, RealType p, RealType q, const Policy& pol)
           {
              const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
              if((p == 0) || (q == 0))
              {
                 //
                 // Can't do a thing if one of p and q is zero:
                 //
                 return policies::raise_evaluation_error<RealType>(function,
                    "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
                    RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
              }
              non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
              tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
              boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
              //
              // Pick an initial guess that we know will give us a probability
              // right around 0.5.
              //
              RealType guess = x - v;
              if(guess < 1)
                 guess = 1;
              std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
                 f, guess, RealType(2), false, tol, max_iter, pol);
              RealType result = ir.first + (ir.second - ir.first) / 2;
              if(max_iter >= policies::get_max_root_iterations<Policy>())
              {
                 policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
                    " or there is no answer to problem.  Current best guess is %1%", result, Policy());
              }
              return result;
           }
  
        }
  
        template <class RealType = double, class Policy = policies::policy<> >
        class non_central_chi_squared_distribution
        {
        public:
           typedef RealType value_type;
           typedef Policy policy_type;
  
           non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
           {
              const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
              RealType r;
              detail::check_df(
                 function,
                 df, &r, Policy());
              detail::check_non_centrality(
                 function,
                 ncp,
                 &r,
                 Policy());
           } // non_central_chi_squared_distribution constructor.
  
           RealType degrees_of_freedom() const
           { // Private data getter function.
              return df;
           }
           RealType non_centrality() const
           { // Private data getter function.
              return ncp;
           }
           static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p)
           {
              const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
              value_type result = detail::find_degrees_of_freedom(
                 static_cast<value_type>(lam),
                 static_cast<value_type>(x),
                 static_cast<value_type>(p),
                 static_cast<value_type>(1-p),
                 forwarding_policy());
              return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 result,
                 function);
           }
           template <class A, class B, class C>
           static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
           {
              const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
              value_type result = detail::find_degrees_of_freedom(
                 static_cast<value_type>(c.dist),
                 static_cast<value_type>(c.param1),
                 static_cast<value_type>(1-c.param2),
                 static_cast<value_type>(c.param2),
                 forwarding_policy());
              return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 result,
                 function);
           }
           static RealType find_non_centrality(RealType v, RealType x, RealType p)
           {
              const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
              value_type result = detail::find_non_centrality(
                 static_cast<value_type>(v),
                 static_cast<value_type>(x),
                 static_cast<value_type>(p),
                 static_cast<value_type>(1-p),
                 forwarding_policy());
              return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 result,
                 function);
           }
           template <class A, class B, class C>
           static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
           {
              const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
              typedef typename policies::evaluation<RealType, Policy>::type value_type;
              typedef typename policies::normalise<
                 Policy,
                 policies::promote_float<false>,
                 policies::promote_double<false>,
                 policies::discrete_quantile<>,
                 policies::assert_undefined<> >::type forwarding_policy;
              value_type result = detail::find_non_centrality(
                 static_cast<value_type>(c.dist),
                 static_cast<value_type>(c.param1),
                 static_cast<value_type>(1-c.param2),
                 static_cast<value_type>(c.param2),
                 forwarding_policy());
              return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                 result,
                 function);
           }
        private:
           // Data member, initialized by constructor.
           RealType df; // degrees of freedom.
           RealType ncp; // non-centrality parameter
        }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
  
        typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
  
        // Non-member functions to give properties of the distribution.
  
        template <class RealType, class Policy>
        inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
        { // Range of permissible values for random variable k.
           using boost::math::tools::max_value;
           return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer?
        }
  
        template <class RealType, class Policy>
        inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
        { // Range of supported values for random variable k.
           // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
           using boost::math::tools::max_value;
           return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
        }
  
        template <class RealType, class Policy>
        inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        { // Mean of poisson distribution = lambda.
           const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy()))
                 return r;
           return k + l;
        } // mean
  
        template <class RealType, class Policy>
        inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        { // mode.
           static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
  
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy()))
                 return (RealType)r;
           return detail::generic_find_mode(dist, 1 + k, function);
        }
  
        template <class RealType, class Policy>
        inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        { // variance.
           const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy()))
                 return r;
           return 2 * (2 * l + k);
        }
  
        // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        // standard_deviation provided by derived accessors.
  
        template <class RealType, class Policy>
        inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        { // skewness = sqrt(l).
           const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy()))
                 return r;
           BOOST_MATH_STD_USING
              return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
        }
  
        template <class RealType, class Policy>
        inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        {
           const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy()))
                 return r;
           return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
        } // kurtosis_excess
  
        template <class RealType, class Policy>
        inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
        {
           return kurtosis_excess(dist) + 3;
        }
  
        template <class RealType, class Policy>
        inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
        { // Probability Density/Mass Function.
           return detail::nccs_pdf(dist, x);
        } // pdf
  
        template <class RealType, class Policy>
        RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
        {
           const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy())
              ||
           !detail::check_positive_x(
              function,
              x,
              &r,
              Policy()))
                 return r;
  
           return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
        } // cdf
  
        template <class RealType, class Policy>
        RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
        { // Complemented Cumulative Distribution Function
           const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
           non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
           RealType x = c.param;
           RealType k = dist.degrees_of_freedom();
           RealType l = dist.non_centrality();
           RealType r;
           if(!detail::check_df(
              function,
              k, &r, Policy())
              ||
           !detail::check_non_centrality(
              function,
              l,
              &r,
              Policy())
              ||
           !detail::check_positive_x(
              function,
              x,
              &r,
              Policy()))
                 return r;
  
           return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
        } // ccdf
  
        template <class RealType, class Policy>
        inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
        { // Quantile (or Percent Point) function.
           return detail::nccs_quantile(dist, p, false);
        } // quantile
  
        template <class RealType, class Policy>
        inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
        { // Quantile (or Percent Point) function.
           return detail::nccs_quantile(c.dist, c.param, true);
        } // quantile complement.
  
     } // namespace math
  } // namespace boost
  
  // This include must be at the end, *after* the accessors
  // for this distribution have been defined, in order to
  // keep compilers that support two-phase lookup happy.
  #include <boost/math/distributions/detail/derived_accessors.hpp>
  
  #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
  
  
  
  