code:
// boost\math\distributions\non_central_t.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_T_HPP
#define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP

#include <boost/math/distributions/fwd.hpp>
#include <boost/math/distributions/non_central_beta.hpp> // for nc beta
#include <boost/math/distributions/normal.hpp> // for normal CDF and quantile
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/detail/generic_quantile.hpp> // quantile

namespace boost
{
   namespace math
   {

      template <class RealType, class Policy>
      class non_central_t_distribution;

      namespace detail{

         template <class T, class Policy>
         T non_central_t2_p(T v, T delta, T x, T y, const Policy& pol, T init_val)
         {
            BOOST_MATH_STD_USING
            //
            // Variables come first:
            //
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
            T errtol = policies::get_epsilon<T, Policy>();
            T d2 = delta * delta / 2;
            //
            // k is the starting point for iteration, and is the
            // maximum of the poisson weighting term, we don't
            // ever allow k == 0 as this can lead to catastrophic
            // cancellation errors later (test case is v = 1621286869049072.3
            // delta = 0.16212868690490723, x = 0.86987415482475994).
            //
            int k = itrunc(d2);
            T pois;
            if(k == 0) k = 1;
            // Starting Poisson weight:
            pois = gamma_p_derivative(T(k+1), d2, pol)
               * tgamma_delta_ratio(T(k + 1), T(0.5f))
               * delta / constants::root_two<T>();
            if(pois == 0)
               return init_val;
            T xterm, beta;
            // Recurrance & starting beta terms:
            beta = x < y
               ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, false, true, &xterm)
               : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, true, true, &xterm);
            xterm *= y / (v / 2 + k);
            T poisf(pois), betaf(beta), xtermf(xterm);
            T sum = init_val;
            if((xterm == 0) && (beta == 0))
               return init_val;

            //
            // Backwards recursion first, this is the stable
            // direction for recursion:
            //
            boost::uintmax_t count = 0;
            T last_term = 0;
            for(int i = k; i >= 0; --i)
            {
               T term = beta * pois;
               sum += term;
               // Don't terminate on first term in case we "fixed" k above:
               if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol)
                  break;
               last_term = term;
               pois *= (i + 0.5f) / d2;
               beta += xterm;
               xterm *= (i) / (x * (v / 2 + i - 1));
               ++count;
            }
            last_term = 0;
            for(int i = k + 1; ; ++i)
            {
               poisf *= d2 / (i + 0.5f);
               xtermf *= (x * (v / 2 + i - 1)) / (i);
               betaf -= xtermf;
               T term = poisf * betaf;
               sum += term;
               if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol))
                  break;
               last_term = term;
               ++count;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "cdf(non_central_t_distribution<%1%>, %1%)",
                     "Series did not converge, closest value was %1%", sum, pol);
               }
            }
            return sum;
         }

         template <class T, class Policy>
         T non_central_t2_q(T v, T delta, T x, T y, const Policy& pol, T init_val)
         {
            BOOST_MATH_STD_USING
            //
            // Variables come first:
            //
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
            T errtol = boost::math::policies::get_epsilon<T, Policy>();
            T d2 = delta * delta / 2;
            //
            // k is the starting point for iteration, and is the
            // maximum of the poisson weighting term, we don't allow
            // k == 0 as this can cause catastrophic cancellation errors
            // (test case is v = 561908036470413.25, delta = 0.056190803647041321,
            // x = 1.6155232703966216):
            //
            int k = itrunc(d2);
            if(k == 0) k = 1;
            // Starting Poisson weight:
            T pois;
            if((k < (int)(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>()))
            {
               //
               // For small k we can optimise this calculation by using
               // a simpler reduced formula:
               //
               pois = exp(-d2);
               pois *= pow(d2, static_cast<T>(k));
               pois /= boost::math::tgamma(T(k + 1 + 0.5), pol);
               pois *= delta / constants::root_two<T>();
            }
            else
            {
               pois = gamma_p_derivative(T(k+1), d2, pol)
                  * tgamma_delta_ratio(T(k + 1), T(0.5f))
                  * delta / constants::root_two<T>();
            }
            if(pois == 0)
               return init_val;
            // Recurance term:
            T xterm;
            T beta;
            // Starting beta term:
            if(k != 0)
            {
               beta = x < y
                  ? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, true, true, &xterm)
                  : detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, false, true, &xterm);

               xterm *= y / (v / 2 + k);
            }
            else
            {
               beta = pow(y, v / 2);
               xterm = beta;
            }
            T poisf(pois), betaf(beta), xtermf(xterm);
            T sum = init_val;
            if((xterm == 0) && (beta == 0))
               return init_val;

            //
            // Fused forward and backwards recursion:
            //
            boost::uintmax_t count = 0;
            T last_term = 0;
            for(int i = k + 1, j = k; ; ++i, --j)
            {
               poisf *= d2 / (i + 0.5f);
               xtermf *= (x * (v / 2 + i - 1)) / (i);
               betaf += xtermf;
               T term = poisf * betaf;

               if(j >= 0)
               {
                  term += beta * pois;
                  pois *= (j + 0.5f) / d2;
                  beta -= xterm;
                  xterm *= (j) / (x * (v / 2 + j - 1));
               }

               sum += term;
               // Don't terminate on first term in case we "fixed" the value of k above:
               if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol)
                  break;
               last_term = term;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "cdf(non_central_t_distribution<%1%>, %1%)",
                     "Series did not converge, closest value was %1%", sum, pol);
               }
               ++count;
            }
            return sum;
         }

         template <class T, class Policy>
         T non_central_t_cdf(T v, T delta, T t, bool invert, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            if ((boost::math::isinf)(v))
            { // Infinite degrees of freedom, so use normal distribution located at delta.
               normal_distribution<T, Policy> n(delta, 1);
               return cdf(n, t);
            }
            //
            // Otherwise, for t < 0 we have to use the reflection formula:
            if(t < 0)
            {
               t = -t;
               delta = -delta;
               invert = !invert;
            }
            if(fabs(delta / (4 * v)) < policies::get_epsilon<T, Policy>())
            {
               // Approximate with a Student's T centred on delta,
               // the crossover point is based on eq 2.6 from
               // "A Comparison of Approximations To Percentiles of the
               // Noncentral t-Distribution".  H. Sahai and M. M. Ojeda,
               // Revista Investigacion Operacional Vol 21, No 2, 2000.
               // Original sources referenced in the above are:
               // "Some Approximations to the Percentage Points of the Noncentral
               // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
               // "Continuous Univariate Distributions".  N.L. Johnson, S. Kotz and
               // N. Balkrishnan. 1995. John Wiley and Sons New York.
               T result = cdf(students_t_distribution<T, Policy>(v), t - delta);
               return invert ? 1 - result : result;
            }
            //
            // x and y are the corresponding random
            // variables for the noncentral beta distribution,
            // with y = 1 - x:
            //
            T x = t * t / (v + t * t);
            T y = v / (v + t * t);
            T d2 = delta * delta;
            T a = 0.5f;
            T b = v / 2;
            T c = a + b + d2 / 2;
            //
            // Crossover point for calculating p or q is the same
            // as for the noncentral beta:
            //
            T cross = 1 - (b / c) * (1 + d2 / (2 * c * c));
            T result;
            if(x < cross)
            {
               //
               // Calculate p:
               //
               if(x != 0)
               {
                  result = non_central_beta_p(a, b, d2, x, y, pol);
                  result = non_central_t2_p(v, delta, x, y, pol, result);
                  result /= 2;
               }
               else
                  result = 0;
               result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
            }
            else
            {
               //
               // Calculate q:
               //
               invert = !invert;
               if(x != 0)
               {
                  result = non_central_beta_q(a, b, d2, x, y, pol);
                  result = non_central_t2_q(v, delta, x, y, pol, result);
                  result /= 2;
               }
               else // x == 0
                  result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta));
            }
            if(invert)
               result = 1 - result;
            return result;
         }

         template <class T, class Policy>
         T non_central_t_quantile(const char* function, T v, T delta, T p, T q, const Policy&)
         {
            BOOST_MATH_STD_USING
     //       static const char* function = "quantile(non_central_t_distribution<%1%>, %1%)";
     // now passed as function
            typedef typename policies::evaluation<T, 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;

               T r;
               if(!detail::check_df_gt0_to_inf(
                  function,
                  v, &r, Policy())
                  ||
               !detail::check_finite(
                  function,
                  delta,
                  &r,
                  Policy())
                  ||
               !detail::check_probability(
                  function,
                  p,
                  &r,
                  Policy()))
                     return r;

            value_type guess = 0;
            if ( ((boost::math::isinf)(v)) || (v > 1 / boost::math::tools::epsilon<T>()) )
            { // Infinite or very large degrees of freedom, so use normal distribution located at delta.
               normal_distribution<T, Policy> n(delta, 1);
               if (p < q)
               {
                 return quantile(n, p);
               }
               else
               {
                 return quantile(complement(n, q));
               }
            }
            else if(v > 3)
            { // Use normal distribution to calculate guess.
               value_type mean = (v > 1 / policies::get_epsilon<T, Policy>()) ? delta : delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f));
               value_type var = (v > 1 / policies::get_epsilon<T, Policy>()) ? value_type(1) : (((delta * delta + 1) * v) / (v - 2) - mean * mean);
               if(p < q)
                  guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
               else
                  guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
            }
            //
            // We *must* get the sign of the initial guess correct,
            // or our root-finder will fail, so double check it now:
            //
            value_type pzero = non_central_t_cdf(
               static_cast<value_type>(v),
               static_cast<value_type>(delta),
               static_cast<value_type>(0),
               !(p < q),
               forwarding_policy());
            int s;
            if(p < q)
               s = boost::math::sign(p - pzero);
            else
               s = boost::math::sign(pzero - q);
            if(s != boost::math::sign(guess))
            {
               guess = static_cast<T>(s);
            }

            value_type result = detail::generic_quantile(
               non_central_t_distribution<value_type, forwarding_policy>(v, delta),
               (p < q ? p : q),
               guess,
               (p >= q),
               function);
            return policies::checked_narrowing_cast<T, forwarding_policy>(
               result,
               function);
         }

         template <class T, class Policy>
         T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val)
         {
            BOOST_MATH_STD_USING
            //
            // Variables come first:
            //
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
            T errtol = boost::math::policies::get_epsilon<T, Policy>();
            T d2 = delta * delta / 2;
            //
            // k is the starting point for iteration, and is the
            // maximum of the poisson weighting term:
            //
            int k = itrunc(d2);
            T pois, xterm;
            if(k == 0)
               k = 1;
            // Starting Poisson weight:
            pois = gamma_p_derivative(T(k+1), d2, pol)
               * tgamma_delta_ratio(T(k + 1), T(0.5f))
               * delta / constants::root_two<T>();
            // Starting beta term:
            xterm = x < y
               ? ibeta_derivative(T(k + 1), n / 2, x, pol)
               : ibeta_derivative(n / 2, T(k + 1), y, pol);
            T poisf(pois), xtermf(xterm);
            T sum = init_val;
            if((pois == 0) || (xterm == 0))
               return init_val;

            //
            // Backwards recursion first, this is the stable
            // direction for recursion:
            //
            boost::uintmax_t count = 0;
            for(int i = k; i >= 0; --i)
            {
               T term = xterm * pois;
               sum += term;
               if(((fabs(term/sum) < errtol) && (i != k)) || (term == 0))
                  break;
               pois *= (i + 0.5f) / d2;
               xterm *= (i) / (x * (n / 2 + i));
               ++count;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "pdf(non_central_t_distribution<%1%>, %1%)",
                     "Series did not converge, closest value was %1%", sum, pol);
               }
            }
            for(int i = k + 1; ; ++i)
            {
               poisf *= d2 / (i + 0.5f);
               xtermf *= (x * (n / 2 + i)) / (i);
               T term = poisf * xtermf;
               sum += term;
               if((fabs(term/sum) < errtol) || (term == 0))
                  break;
               ++count;
               if(count > max_iter)
               {
                  return policies::raise_evaluation_error(
                     "pdf(non_central_t_distribution<%1%>, %1%)",
                     "Series did not converge, closest value was %1%", sum, pol);
               }
            }
            return sum;
         }

         template <class T, class Policy>
         T non_central_t_pdf(T n, T delta, T t, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            if ((boost::math::isinf)(n))
            { // Infinite degrees of freedom, so use normal distribution located at delta.
               normal_distribution<T, Policy> norm(delta, 1);
               return pdf(norm, t);
            }
            //
            // Otherwise, for t < 0 we have to use the reflection formula:
            if(t < 0)
            {
               t = -t;
               delta = -delta;
            }
            if(t == 0)
            {
               //
               // Handle this as a special case, using the formula
               // from Weisstein, Eric W.
               // "Noncentral Student's t-Distribution."
               // From MathWorld--A Wolfram Web Resource.
               // http://mathworld.wolfram.com/NoncentralStudentst-Distribution.html
               //
               // The formula is simplified thanks to the relation
               // 1F1(a,b,0) = 1.
               //
               return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f))
                  * sqrt(n / constants::pi<T>())
                  * exp(-delta * delta / 2) / 2;
            }
            if(fabs(delta / (4 * n)) < policies::get_epsilon<T, Policy>())
            {
               // Approximate with a Student's T centred on delta,
               // the crossover point is based on eq 2.6 from
               // "A Comparison of Approximations To Percentiles of the
               // Noncentral t-Distribution".  H. Sahai and M. M. Ojeda,
               // Revista Investigacion Operacional Vol 21, No 2, 2000.
               // Original sources referenced in the above are:
               // "Some Approximations to the Percentage Points of the Noncentral
               // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
               // "Continuous Univariate Distributions".  N.L. Johnson, S. Kotz and
               // N. Balkrishnan. 1995. John Wiley and Sons New York.
               return pdf(students_t_distribution<T, Policy>(n), t - delta);
            }
            //
            // x and y are the corresponding random
            // variables for the noncentral beta distribution,
            // with y = 1 - x:
            //
            T x = t * t / (n + t * t);
            T y = n / (n + t * t);
            T a = 0.5f;
            T b = n / 2;
            T d2 = delta * delta;
            //
            // Calculate pdf:
            //
            T dt = n * t / (n * n + 2 * n * t * t + t * t * t * t);
            T result = non_central_beta_pdf(a, b, d2, x, y, pol);
            T tol = tools::epsilon<T>() * result * 500;
            result = non_central_t2_pdf(n, delta, x, y, pol, result);
            if(result <= tol)
               result = 0;
            result *= dt;
            return result;
         }

         template <class T, class Policy>
         T mean(T v, T delta, const Policy& pol)
         {
            if ((boost::math::isinf)(v))
            {
               return delta;
            }
            BOOST_MATH_STD_USING
            if (v > 1 / boost::math::tools::epsilon<T>() )
            {
              //normal_distribution<T, Policy> n(delta, 1);
              //return boost::math::mean(n);
              return delta;
            }
            else
            {
             return delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol);
            }
            // Other moments use mean so using normal distribution is propagated.
         }

         template <class T, class Policy>
         T variance(T v, T delta, const Policy& pol)
         {
            if ((boost::math::isinf)(v))
            {
               return 1;
            }
            if (delta == 0)
            {  // == Student's t
              return v / (v - 2);
            }
            T result = ((delta * delta + 1) * v) / (v - 2);
            T m = mean(v, delta, pol);
            result -= m * m;
            return result;
         }

         template <class T, class Policy>
         T skewness(T v, T delta, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            if ((boost::math::isinf)(v))
            {
               return 0;
            }
            if(delta == 0)
            { // == Student's t
              return 0;
            }
            T mean = boost::math::detail::mean(v, delta, pol);
            T l2 = delta * delta;
            T var = ((l2 + 1) * v) / (v - 2) - mean * mean;
            T result = -2 * var;
            result += v * (l2 + 2 * v - 3) / ((v - 3) * (v - 2));
            result *= mean;
            result /= pow(var, T(1.5f));
            return result;
         }

         template <class T, class Policy>
         T kurtosis_excess(T v, T delta, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            if ((boost::math::isinf)(v))
            {
               return 3;
            }
            if (delta == 0)
            { // == Student's t
              return 3;
            }
            T mean = boost::math::detail::mean(v, delta, pol);
            T l2 = delta * delta;
            T var = ((l2 + 1) * v) / (v - 2) - mean * mean;
            T result = -3 * var;
            result += v * (l2 * (v + 1) + 3 * (3 * v - 5)) / ((v - 3) * (v - 2));
            result *= -mean * mean;
            result += v * v * (l2 * l2 + 6 * l2 + 3) / ((v - 4) * (v - 2));
            result /= var * var;
            return result;
         }

#if 0
         //
         // This code is disabled, since there can be multiple answers to the
         // question, and it's not clear how to find the "right" one.
         //
         template <class RealType, class Policy>
         struct t_degrees_of_freedom_finder
         {
            t_degrees_of_freedom_finder(
               RealType delta_, RealType x_, RealType p_, bool c)
               : delta(delta_), x(x_), p(p_), comp(c) {}

            RealType operator()(const RealType& v)
            {
               non_central_t_distribution<RealType, Policy> d(v, delta);
               return comp ?
                  p - cdf(complement(d, x))
                  : cdf(d, x) - p;
            }
         private:
            RealType delta;
            RealType x;
            RealType p;
            bool comp;
         };

         template <class RealType, class Policy>
         inline RealType find_t_degrees_of_freedom(
            RealType delta, RealType x, RealType p, RealType q, const Policy& pol)
         {
            const char* function = "non_central_t<%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());
            }
            t_degrees_of_freedom_finder<RealType, Policy> f(delta, 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:
            //
            RealType guess = 200;
            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>())
            {
               return 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 t_non_centrality_finder
         {
            t_non_centrality_finder(
               RealType v_, RealType x_, RealType p_, bool c)
               : v(v_), x(x_), p(p_), comp(c) {}

            RealType operator()(const RealType& delta)
            {
               non_central_t_distribution<RealType, Policy> d(v, delta);
               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_t_non_centrality(
            RealType v, RealType x, RealType p, RealType q, const Policy& pol)
         {
            const char* function = "non_central_t<%1%>::find_t_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());
            }
            t_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 is the right side of
            // zero:
            //
            RealType guess;
            if(f(0) < 0)
               guess = 1;
            else
               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>())
            {
               return 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;
         }
#endif
      } // namespace detail ======================================================================

      template <class RealType = double, class Policy = policies::policy<> >
      class non_central_t_distribution
      {
      public:
         typedef RealType value_type;
         typedef Policy policy_type;

         non_central_t_distribution(RealType v_, RealType lambda) : v(v_), ncp(lambda)
         {
            const char* function = "boost::math::non_central_t_distribution<%1%>::non_central_t_distribution(%1%,%1%)";
            RealType r;
            detail::check_df_gt0_to_inf(
               function,
               v, &r, Policy());
            detail::check_finite(
               function,
               lambda,
               &r,
               Policy());
         } // non_central_t_distribution constructor.

         RealType degrees_of_freedom() const
         { // Private data getter function.
            return v;
         }
         RealType non_centrality() const
         { // Private data getter function.
            return ncp;
         }
#if 0
         //
         // This code is disabled, since there can be multiple answers to the
         // question, and it's not clear how to find the "right" one.
         //
         static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p)
         {
            const char* function = "non_central_t<%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_t_degrees_of_freedom(
               static_cast<value_type>(delta),
               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_t<%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_t_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_t<%1%>::find_t_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_t_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_t<%1%>::find_t_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_t_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);
         }
#endif
      private:
         // Data member, initialized by constructor.
         RealType v;   // degrees of freedom
         RealType ncp; // non-centrality parameter
      }; // template <class RealType, class Policy> class non_central_t_distribution

      typedef non_central_t_distribution<double> non_central_t; // 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_t_distribution<RealType, Policy>& /* dist */)
      { // Range of permissible values for random variable k.
         using boost::math::tools::max_value;
         return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
      }

      template <class RealType, class Policy>
      inline const std::pair<RealType, RealType> support(const non_central_t_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>(-max_value<RealType>(), max_value<RealType>());
      }

      template <class RealType, class Policy>
      inline RealType mode(const non_central_t_distribution<RealType, Policy>& dist)
      { // mode.
         static const char* function = "mode(non_central_t_distribution<%1%> const&)";
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy()))
               return (RealType)r;

         BOOST_MATH_STD_USING

         RealType m = v < 3 ? 0 : detail::mean(v, l, Policy());
         RealType var = v < 4 ? 1 : detail::variance(v, l, Policy());

         return detail::generic_find_mode(
            dist,
            m,
            function,
            sqrt(var));
      }

      template <class RealType, class Policy>
      inline RealType mean(const non_central_t_distribution<RealType, Policy>& dist)
      {
         BOOST_MATH_STD_USING
         const char* function = "mean(const non_central_t_distribution<%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;
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy()))
               return (RealType)r;
         if(v <= 1)
            return policies::raise_domain_error<RealType>(
               function,
               "The non-central t distribution has no defined mean for degrees of freedom <= 1: got v=%1%.", v, Policy());
         // return l * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, RealType(0.5f));
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::mean(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);

      } // mean

      template <class RealType, class Policy>
      inline RealType variance(const non_central_t_distribution<RealType, Policy>& dist)
      { // variance.
         const char* function = "variance(const non_central_t_distribution<%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;
         BOOST_MATH_STD_USING
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy()))
               return (RealType)r;
         if(v <= 2)
            return policies::raise_domain_error<RealType>(
               function,
               "The non-central t distribution has no defined variance for degrees of freedom <= 2: got v=%1%.", v, Policy());
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::variance(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
      }

      // RealType standard_deviation(const non_central_t_distribution<RealType, Policy>& dist)
      // standard_deviation provided by derived accessors.

      template <class RealType, class Policy>
      inline RealType skewness(const non_central_t_distribution<RealType, Policy>& dist)
      { // skewness = sqrt(l).
         const char* function = "skewness(const non_central_t_distribution<%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;
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy()))
               return (RealType)r;
         if(v <= 3)
            return policies::raise_domain_error<RealType>(
               function,
               "The non-central t distribution has no defined skewness for degrees of freedom <= 3: got v=%1%.", v, Policy());;
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::skewness(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
      }

      template <class RealType, class Policy>
      inline RealType kurtosis_excess(const non_central_t_distribution<RealType, Policy>& dist)
      {
         const char* function = "kurtosis_excess(const non_central_t_distribution<%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;
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy()))
               return (RealType)r;
         if(v <= 4)
            return policies::raise_domain_error<RealType>(
               function,
               "The non-central t distribution has no defined kurtosis for degrees of freedom <= 4: got v=%1%.", v, Policy());;
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::kurtosis_excess(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
      } // kurtosis_excess

      template <class RealType, class Policy>
      inline RealType kurtosis(const non_central_t_distribution<RealType, Policy>& dist)
      {
         return kurtosis_excess(dist) + 3;
      }

      template <class RealType, class Policy>
      inline RealType pdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& t)
      { // Probability Density/Mass Function.
         const char* function = "pdf(non_central_t_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;

         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy())
            ||
         !detail::check_x(
            function,
            t,
            &r,
            Policy()))
               return (RealType)r;
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::non_central_t_pdf(static_cast<value_type>(v),
               static_cast<value_type>(l),
               static_cast<value_type>(t),
               Policy()),
            function);
      } // pdf

      template <class RealType, class Policy>
      RealType cdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& x)
      {
         const char* function = "boost::math::cdf(non_central_t_distribution<%1%>&, %1%)";
//   was const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%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;

         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy())
            ||
         !detail::check_x(
            function,
            x,
            &r,
            Policy()))
               return (RealType)r;
          if ((boost::math::isinf)(v))
          { // Infinite degrees of freedom, so use normal distribution located at delta.
             normal_distribution<RealType, Policy> n(l, 1);
             cdf(n, x);
              //return cdf(normal_distribution<RealType, Policy>(l, 1), x);
          }

         if(l == 0)
         { // NO non-centrality, so use Student's t instead.
            return cdf(students_t_distribution<RealType, Policy>(v), x);
         }
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::non_central_t_cdf(
               static_cast<value_type>(v),
               static_cast<value_type>(l),
               static_cast<value_type>(x),
               false, Policy()),
            function);
      } // cdf

      template <class RealType, class Policy>
      RealType cdf(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c)
      { // Complemented Cumulative Distribution Function
  // was       const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)";
         const char* function = "boost::math::cdf(const complement(non_central_t_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;

         non_central_t_distribution<RealType, Policy> const& dist = c.dist;
         RealType x = c.param;
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality(); // aka delta
         RealType r;
         if(!detail::check_df_gt0_to_inf(
            function,
            v, &r, Policy())
            ||
         !detail::check_finite(
            function,
            l,
            &r,
            Policy())
            ||
         !detail::check_x(
            function,
            x,
            &r,
            Policy()))
               return (RealType)r;

         if ((boost::math::isinf)(v))
         { // Infinite degrees of freedom, so use normal distribution located at delta.
             normal_distribution<RealType, Policy> n(l, 1);
             return cdf(complement(n, x));
         }
         if(l == 0)
         { // zero non-centrality so use Student's t distribution.
            return cdf(complement(students_t_distribution<RealType, Policy>(v), x));
         }
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
            detail::non_central_t_cdf(
               static_cast<value_type>(v),
               static_cast<value_type>(l),
               static_cast<value_type>(x),
               true, Policy()),
            function);
      } // ccdf

      template <class RealType, class Policy>
      inline RealType quantile(const non_central_t_distribution<RealType, Policy>& dist, const RealType& p)
      { // Quantile (or Percent Point) function.
         static const char* function = "quantile(const non_central_t_distribution<%1%>, %1%)";
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         return detail::non_central_t_quantile(function, v, l, p, RealType(1-p), Policy());
      } // quantile

      template <class RealType, class Policy>
      inline RealType quantile(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c)
      { // Quantile (or Percent Point) function.
         static const char* function = "quantile(const complement(non_central_t_distribution<%1%>, %1%))";
         non_central_t_distribution<RealType, Policy> const& dist = c.dist;
         RealType q = c.param;
         RealType v = dist.degrees_of_freedom();
         RealType l = dist.non_centrality();
         return detail::non_central_t_quantile(function, v, l, RealType(1-q), q, Policy());
      } // 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_T_HPP


code:
//  Boost utility.hpp header file  -------------------------------------------//

//  Copyright 1999-2003 Aleksey Gurtovoy.  Use, modification, and distribution are
//  subject to the Boost Software License, Version 1.0.  (See accompanying file
//  LICENSE_1_0.txt or a copy at <http://www.boost.org/LICENSE_1_0.txt>.)

//  See <http://www.boost.org/libs/utility/> for the library's home page.

#ifndef BOOST_UTILITY_HPP
#define BOOST_UTILITY_HPP

#include <boost/utility/addressof.hpp>
#include <boost/utility/base_from_member.hpp>
#include <boost/utility/binary.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/utility/identity_type.hpp>
#include <boost/checked_delete.hpp>
#include <boost/next_prior.hpp>
#include <boost/noncopyable.hpp>

#endif  // BOOST_UTILITY_HPP

license: mit


code:
//  (C) Copyright John Maddock 2006.
//  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_LEGENDRE_HPP
#define BOOST_MATH_SPECIAL_LEGENDRE_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/tools/config.hpp>

namespace boost{
namespace math{

// Recurrance relation for legendre P and Q polynomials:
template <class T1, class T2, class T3>
inline typename tools::promote_args<T1, T2, T3>::type
   legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1)
{
   typedef typename tools::promote_args<T1, T2, T3>::type result_type;
   return ((2 * l + 1) * result_type(x) * result_type(Pl) - l * result_type(Plm1)) / (l + 1);
}

namespace detail{

// Implement Legendre P and Q polynomials via recurrance:
template <class T, class Policy>
T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)
{
   static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)";
   // Error handling:
   if((x < -1) || (x > 1))
      return policies::raise_domain_error<T>(
         function,
         "The Legendre Polynomial is defined for"
         " -1 <= x <= 1, but got x = %1%.", x, pol);

   T p0, p1;
   if(second)
   {
      // A solution of the second kind (Q):
      p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;
      p1 = x * p0 - 1;
   }
   else
   {
      // A solution of the first kind (P):
      p0 = 1;
      p1 = x;
   }
   if(l == 0)
      return p0;

   unsigned n = 1;

   while(n < l)
   {
      std::swap(p0, p1);
      p1 = boost::math::legendre_next(n, x, p0, p1);
      ++n;
   }
   return p1;
}

} // namespace detail

template <class T, class Policy>
inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
   legendre_p(int l, T x, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   static const char* function = "boost::math::legendre_p<%1%>(unsigned, %1%)";
   if(l < 0)
      return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(-l-1, static_cast<value_type>(x), pol, false), function);
   return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function);
}

template <class T>
inline typename tools::promote_args<T>::type
   legendre_p(int l, T x)
{
   return boost::math::legendre_p(l, x, policies::policy<>());
}

template <class T, class Policy>
inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
   legendre_q(unsigned l, T x, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, true), "boost::math::legendre_q<%1%>(unsigned, %1%)");
}

template <class T>
inline typename tools::promote_args<T>::type
   legendre_q(unsigned l, T x)
{
   return boost::math::legendre_q(l, x, policies::policy<>());
}

// Recurrence for associated polynomials:
template <class T1, class T2, class T3>
inline typename tools::promote_args<T1, T2, T3>::type
   legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1)
{
   typedef typename tools::promote_args<T1, T2, T3>::type result_type;
   return ((2 * l + 1) * result_type(x) * result_type(Pl) - (l + m) * result_type(Plm1)) / (l + 1 - m);
}

namespace detail{
// Legendre P associated polynomial:
template <class T, class Policy>
T legendre_p_imp(int l, int m, T x, T sin_theta_power, const Policy& pol)
{
   // Error handling:
   if((x < -1) || (x > 1))
      return policies::raise_domain_error<T>(
      "boost::math::legendre_p<%1%>(int, int, %1%)",
         "The associated Legendre Polynomial is defined for"
         " -1 <= x <= 1, but got x = %1%.", x, pol);
   // Handle negative arguments first:
   if(l < 0)
      return legendre_p_imp(-l-1, m, x, sin_theta_power, pol);
   if(m < 0)
   {
      int sign = (m&1) ? -1 : 1;
      return sign * boost::math::tgamma_ratio(static_cast<T>(l+m+1), static_cast<T>(l+1-m), pol) * legendre_p_imp(l, -m, x, sin_theta_power, pol);
   }
   // Special cases:
   if(m > l)
      return 0;
   if(m == 0)
      return boost::math::legendre_p(l, x, pol);

   T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power;

   if(m&1)
      p0 *= -1;
   if(m == l)
      return p0;

   T p1 = x * (2 * m + 1) * p0;

   int n = m + 1;

   while(n < l)
   {
      std::swap(p0, p1);
      p1 = boost::math::legendre_next(n, m, x, p0, p1);
      ++n;
   }
   return p1;
}

template <class T, class Policy>
inline T legendre_p_imp(int l, int m, T x, const Policy& pol)
{
   BOOST_MATH_STD_USING
   // TODO: we really could use that mythical "pow1p" function here:
   return legendre_p_imp(l, m, x, static_cast<T>(pow(1 - x*x, T(abs(m))/2)), pol);
}

}

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   legendre_p(int l, int m, T x, const Policy& pol)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_imp(l, m, static_cast<value_type>(x), pol), "bost::math::legendre_p<%1%>(int, int, %1%)");
}

template <class T>
inline typename tools::promote_args<T>::type
   legendre_p(int l, int m, T x)
{
   return boost::math::legendre_p(l, m, x, policies::policy<>());
}

} // namespace math
} // namespace boost

#endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP


code:
//---------------------------------------------------------------------------//
// Copyright (c) 2013-2014 Kyle Lutz <kyle.r.lutz@gmail.com>
//
// Distributed under 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
//
// See http://boostorg.github.com/compute for more information.
//---------------------------------------------------------------------------//

#ifndef BOOST_COMPUTE_ITERATOR_DISCARD_ITERATOR_HPP
#define BOOST_COMPUTE_ITERATOR_DISCARD_ITERATOR_HPP

#include <string>
#include <cstddef>
#include <iterator>

#include <boost/config.hpp>
#include <boost/iterator/iterator_facade.hpp>

#include <boost/compute/detail/meta_kernel.hpp>
#include <boost/compute/type_traits/is_device_iterator.hpp>

namespace boost {
namespace compute {

// forward declaration for discard_iterator
class discard_iterator;

namespace detail {

// helper class which defines the iterator_facade super-class
// type for discard_iterator
struct discard_iterator_base
{
    typedef ::boost::iterator_facade<
        ::boost::compute::discard_iterator,
        void,
        ::std::random_access_iterator_tag,
        void *
    > type;
};

template<class IndexExpr>
struct discard_iterator_index_expr
{
    typedef void result_type;

    discard_iterator_index_expr(const IndexExpr &expr)
        : m_expr(expr)
    {
    }

    IndexExpr m_expr;
};

template<class IndexExpr>
inline meta_kernel& operator<<(meta_kernel &kernel,
                               const discard_iterator_index_expr<IndexExpr> &expr)
{
    (void) expr;

    return kernel;
}

} // end detail namespace

/// \class discard_iterator
/// \brief An iterator which discards all values written to it.
///
/// \see make_discard_iterator(), constant_iterator
class discard_iterator : public detail::discard_iterator_base::type
{
public:
    typedef detail::discard_iterator_base::type super_type;
    typedef super_type::reference reference;
    typedef super_type::difference_type difference_type;

    discard_iterator(size_t index = 0)
        : m_index(index)
    {
    }

    discard_iterator(const discard_iterator &other)
        : m_index(other.m_index)
    {
    }

    discard_iterator& operator=(const discard_iterator &other)
    {
        if(this != &other){
            m_index = other.m_index;
        }

        return *this;
    }

    ~discard_iterator()
    {
    }

    /// \internal_
    template<class Expr>
    detail::discard_iterator_index_expr<Expr>
    operator[](const Expr &expr) const
    {
        return detail::discard_iterator_index_expr<Expr>(expr);
    }

private:
    friend class ::boost::iterator_core_access;

    /// \internal_
    reference dereference() const
    {
        return 0;
    }

    /// \internal_
    bool equal(const discard_iterator &other) const
    {
        return m_index == other.m_index;
    }

    /// \internal_
    void increment()
    {
        m_index++;
    }

    /// \internal_
    void decrement()
    {
        m_index--;
    }

    /// \internal_
    void advance(difference_type n)
    {
        m_index = static_cast<size_t>(static_cast<difference_type>(m_index) + n);
    }

    /// \internal_
    difference_type distance_to(const discard_iterator &other) const
    {
        return static_cast<difference_type>(other.m_index - m_index);
    }

private:
    size_t m_index;
};

/// Returns a new discard_iterator with \p index.
///
/// \param index the index of the iterator
///
/// \return a \c discard_iterator at \p index
inline discard_iterator make_discard_iterator(size_t index = 0)
{
    return discard_iterator(index);
}

/// internal_ (is_device_iterator specialization for discard_iterator)
template<>
struct is_device_iterator<discard_iterator> : boost::true_type {};

} // end compute namespace
} // end boost namespace

#endif // BOOST_COMPUTE_ITERATOR_DISCARD_ITERATOR_HPP


code:
#ifndef BOOST_MPL_BITXOR_HPP_INCLUDED
#define BOOST_MPL_BITXOR_HPP_INCLUDED

// Copyright Aleksey Gurtovoy 2000-2004
// Copyright Jaap Suter 2003
//
// Distributed under 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)
//
// See http://www.boost.org/libs/mpl for documentation.

// $Id$
// $Date$
// $Revision$

#define AUX778076_OP_NAME   bitxor_
#define AUX778076_OP_PREFIX bitxor
#define AUX778076_OP_TOKEN  ^
#include <boost/mpl/aux_/arithmetic_op.hpp>

#endif // BOOST_MPL_BITXOR_HPP_INCLUDED


code:
#ifndef DATE_TIME_TIME_ITERATOR_HPP___
#define DATE_TIME_TIME_ITERATOR_HPP___

/* Copyright (c) 2002,2003 CrystalClear Software, Inc.
 * Use, modification and distribution is subject to the
 * Boost Software License, Version 1.0. (See accompanying
 * file LICENSE_1_0.txt or http://www.boost.org/LICENSE_1_0.txt)
 * Author: Jeff Garland, Bart Garst
 * $Date$
 */

namespace boost {
namespace date_time {

  //! Simple time iterator skeleton class
  template<class time_type>
  class time_itr {
  public:
    typedef typename time_type::time_duration_type time_duration_type;
    time_itr(time_type t, time_duration_type d) : current_(t), offset_(d) {}
    time_itr& operator++()
    {
      current_ = current_ + offset_;
      return *this;
    }
    time_itr& operator--()
    {
      current_ = current_ - offset_;
      return *this;
    }
    time_type operator*() {return current_;}
    time_type* operator->() {return &current_;}
    bool operator<  (const time_type& t) {return current_ < t;}
    bool operator<= (const time_type& t) {return current_ <= t;}
    bool operator!=  (const time_type& t) {return current_ != t;}
    bool operator== (const time_type& t) {return current_ == t;}
    bool operator>  (const time_type& t) {return current_ > t;}
    bool operator>= (const time_type& t) {return current_ >= t;}

  private:
    time_type current_;
    time_duration_type offset_;
  };

} }//namespace date_time

#endif


code:
//  (C) Copyright R.W. Grosse-Kunstleve 2002.
//  Distributed under 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_PYTHON_OBJECT_PICKLE_SUPPORT_RWGK20020603_HPP
# define BOOST_PYTHON_OBJECT_PICKLE_SUPPORT_RWGK20020603_HPP

# include <boost/python/detail/prefix.hpp>

namespace boost { namespace python {

namespace api
{
  class object;
}
using api::object;
class tuple;

BOOST_PYTHON_DECL object const& make_instance_reduce_function();

struct pickle_suite;

namespace error_messages {

  template <class T>
  struct missing_pickle_suite_function_or_incorrect_signature {};

  inline void must_be_derived_from_pickle_suite(pickle_suite const&) {}
}

namespace detail { struct pickle_suite_registration; }

struct pickle_suite
{
  private:
    struct inaccessible {};
    friend struct detail::pickle_suite_registration;
  public:
    static inaccessible* getinitargs() { return 0; }
    static inaccessible* getstate() { return 0; }
    static inaccessible* setstate() { return 0; }
    static bool getstate_manages_dict() { return false; }
};

namespace detail {

  struct pickle_suite_registration
  {
    typedef pickle_suite::inaccessible inaccessible;

    template <class Class_, class Tgetinitargs>
    static
    void
    register_(
      Class_& cl,
      tuple (*getinitargs_fn)(Tgetinitargs),
      inaccessible* (* /*getstate_fn*/)(),
      inaccessible* (* /*setstate_fn*/)(),
      bool)
    {
      cl.enable_pickling_(false);
      cl.def("__getinitargs__", getinitargs_fn);
    }

    template <class Class_,
              class Rgetstate, class Tgetstate,
              class Tsetstate, class Ttuple>
    static
    void
    register_(
      Class_& cl,
      inaccessible* (* /*getinitargs_fn*/)(),
      Rgetstate (*getstate_fn)(Tgetstate),
      void (*setstate_fn)(Tsetstate, Ttuple),
      bool getstate_manages_dict)
    {
      cl.enable_pickling_(getstate_manages_dict);
      cl.def("__getstate__", getstate_fn);
      cl.def("__setstate__", setstate_fn);
    }

    template <class Class_,
              class Tgetinitargs,
              class Rgetstate, class Tgetstate,
              class Tsetstate, class Ttuple>
    static
    void
    register_(
      Class_& cl,
      tuple (*getinitargs_fn)(Tgetinitargs),
      Rgetstate (*getstate_fn)(Tgetstate),
      void (*setstate_fn)(Tsetstate, Ttuple),
      bool getstate_manages_dict)
    {
      cl.enable_pickling_(getstate_manages_dict);
      cl.def("__getinitargs__", getinitargs_fn);
      cl.def("__getstate__", getstate_fn);
      cl.def("__setstate__", setstate_fn);
    }

    template <class Class_>
    static
    void
    register_(
      Class_&,
      ...)
    {
      typedef typename
        error_messages::missing_pickle_suite_function_or_incorrect_signature<
          Class_>::error_type error_type BOOST_ATTRIBUTE_UNUSED;
    }
  };

  template <typename PickleSuiteType>
  struct pickle_suite_finalize
  : PickleSuiteType,
    pickle_suite_registration
  {};

} // namespace detail

}} // namespace boost::python

#endif // BOOST_PYTHON_OBJECT_PICKLE_SUPPORT_RWGK20020603_HPP


code:
#ifndef BOOST_THREAD_DETAIL_THREAD_GROUP_HPP
#define BOOST_THREAD_DETAIL_THREAD_GROUP_HPP
// Distributed under 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)
// (C) Copyright 2007-9 Anthony Williams

#include <list>
#include <boost/thread/csbl/memory/unique_ptr.hpp>
#include <boost/thread/shared_mutex.hpp>
#include <boost/thread/mutex.hpp>
#include <boost/thread/lock_guard.hpp>

#include <boost/config/abi_prefix.hpp>

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable:4251)
#endif

namespace boost
{
    class thread_group
    {
    private:
        thread_group(thread_group const&);
        thread_group& operator=(thread_group const&);
    public:
        thread_group() {}
        ~thread_group()
        {
            for(std::list<thread*>::iterator it=threads.begin(),end=threads.end();
                it!=end;
                ++it)
            {
                delete *it;
            }
        }

        bool is_this_thread_in()
        {
            thread::id id = this_thread::get_id();
            boost::shared_lock<shared_mutex> guard(m);
            for(std::list<thread*>::iterator it=threads.begin(),end=threads.end();
                it!=end;
                ++it)
            {
              if ((*it)->get_id() == id)
                return true;
            }
            return false;
        }

        bool is_thread_in(thread* thrd)
        {
          if(thrd)
          {
            thread::id id = thrd->get_id();
            boost::shared_lock<shared_mutex> guard(m);
            for(std::list<thread*>::iterator it=threads.begin(),end=threads.end();
                it!=end;
                ++it)
            {
              if ((*it)->get_id() == id)
                return true;
            }
            return false;
          }
          else
          {
            return false;
          }
        }

        template<typename F>
        thread* create_thread(F threadfunc)
        {
            boost::lock_guard<shared_mutex> guard(m);
            boost::csbl::unique_ptr<thread> new_thread(new thread(threadfunc));
            threads.push_back(new_thread.get());
            return new_thread.release();
        }

        void add_thread(thread* thrd)
        {
            if(thrd)
            {
                BOOST_THREAD_ASSERT_PRECONDITION( ! is_thread_in(thrd) ,
                    thread_resource_error(static_cast<int>(system::errc::resource_deadlock_would_occur), "boost::thread_group: trying to add a duplicated thread")
                );

                boost::lock_guard<shared_mutex> guard(m);
                threads.push_back(thrd);
            }
        }

        void remove_thread(thread* thrd)
        {
            boost::lock_guard<shared_mutex> guard(m);
            std::list<thread*>::iterator const it=std::find(threads.begin(),threads.end(),thrd);
            if(it!=threads.end())
            {
                threads.erase(it);
            }
        }

        void join_all()
        {
            BOOST_THREAD_ASSERT_PRECONDITION( ! is_this_thread_in() ,
                thread_resource_error(static_cast<int>(system::errc::resource_deadlock_would_occur), "boost::thread_group: trying joining itself")
            );
            boost::shared_lock<shared_mutex> guard(m);

            for(std::list<thread*>::iterator it=threads.begin(),end=threads.end();
                it!=end;
                ++it)
            {
              if ((*it)->joinable())
                (*it)->join();
            }
        }

#if defined BOOST_THREAD_PROVIDES_INTERRUPTIONS
        void interrupt_all()
        {
            boost::shared_lock<shared_mutex> guard(m);

            for(std::list<thread*>::iterator it=threads.begin(),end=threads.end();
                it!=end;
                ++it)
            {
                (*it)->interrupt();
            }
        }
#endif

        size_t size() const
        {
            boost::shared_lock<shared_mutex> guard(m);
            return threads.size();
        }

    private:
        std::list<thread*> threads;
        mutable shared_mutex m;
    };
}

#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

#include <boost/config/abi_suffix.hpp>

#endif


code:
//  (C) Copyright Gennadiy Rozental 2001.
//  Distributed under 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)

//  See http://www.boost.org/libs/test for the library home page.
//
/// @file
/// @brief Defines unit test log formatter interface
///
/// You can define a class with implements this interface and use an instance of it
/// as a Unit Test Framework log formatter
// ***************************************************************************

#ifndef BOOST_TEST_UNIT_TEST_LOG_FORMATTER_HPP_071894GER
#define BOOST_TEST_UNIT_TEST_LOG_FORMATTER_HPP_071894GER

// Boost.Test
#include <boost/test/detail/global_typedef.hpp>
#include <boost/test/detail/log_level.hpp>
#include <boost/test/detail/fwd_decl.hpp>

// STL
#include <iosfwd>
#include <string> // for std::string
#include <iostream>

#include <boost/test/detail/suppress_warnings.hpp>

//____________________________________________________________________________//

namespace boost {
namespace unit_test {

// ************************************************************************** //
/// Collection of log entry attributes
// ************************************************************************** //

struct BOOST_TEST_DECL log_entry_data {
    log_entry_data()
    {
        m_file_name.reserve( 200 );
    }

    std::string     m_file_name; ///< log entry file name
    std::size_t     m_line_num;  ///< log entry line number
    log_level       m_level;     ///< log entry level

    void clear()
    {
        m_file_name.erase();
        m_line_num      = 0;
        m_level     = log_nothing;
    }
};

// ************************************************************************** //
/// Collection of log checkpoint attributes
// ************************************************************************** //

struct BOOST_TEST_DECL log_checkpoint_data
{
    const_string    m_file_name; ///< log checkpoint file name
    std::size_t     m_line_num;  ///< log checkpoint file name
    std::string     m_message;   ///< log checkpoint message

    void clear()
    {
        m_file_name.clear();
        m_line_num  = 0;
        m_message   = std::string();
    }
};

// ************************************************************************** //
/// @brief Abstract Unit Test Framework log formatter interface
///
/// During the test module execution Unit Test Framework can report messages about success
/// or failure of assertions, which test suites are being run and more (specifically which
/// messages are reported depends on log level threshold selected by the user).
///
/// All these messages constitute Unit Test Framework log. There are many ways (formats) to present
/// these messages to the user.
///
/// Boost.Test comes with three formats:
/// - Compiler-like log format: intended for human consumption/diagnostic
/// - XML based log format:  intended for processing by automated regression test systems.
/// - JUNIT based log format:  intended for processing by automated regression test systems.
///
/// If you want to produce some other format you need to implement class with specific interface and use
/// method @c unit_test_log_t::set_formatter during a test module initialization to set an active formatter.
/// The class unit_test_log_formatter defines this interface.
///
/// This interface requires you to format all possible messages being produced in the log.
/// These includes error messages about failed assertions, messages about caught exceptions and
/// information messages about test units being started/ended. All the methods in this interface takes
/// a reference to standard stream as a first argument. This is where final messages needs to be directed
/// to. Also you are given all the information necessary to produce a message.
///
/// @par Since Boost 1.62:
/// - Each formatter may indicate the default output stream. This is convenient for instance for streams intended
///   for automated processing that indicate a file. See @c get_default_stream_description for more details.
/// - Each formatter may manage its own log level through the getter/setter @c get_log_level and @c set_log_level .
///
/// @see
/// - boost::unit_test::test_observer for an indication of the calls of the test observer interface
class BOOST_TEST_DECL unit_test_log_formatter {
public:
    /// Types of log entries (messages written into a log)
    enum log_entry_types { BOOST_UTL_ET_INFO,       ///< Information message from the framework
                           BOOST_UTL_ET_MESSAGE,    ///< Information message from the user
                           BOOST_UTL_ET_WARNING,    ///< Warning (non error) condition notification message
                           BOOST_UTL_ET_ERROR,      ///< Non fatal error notification message
                           BOOST_UTL_ET_FATAL_ERROR ///< Fatal error notification message
    };

    //! Constructor
    unit_test_log_formatter()
        : m_log_level(log_all_errors)
    {}

    // Destructor
    virtual             ~unit_test_log_formatter() {}

    // @name Test start/finish

    /// Invoked at the beginning of test module execution
    ///
    /// @param[in] os   output stream to write a messages to
    /// @param[in] test_cases_amount total test case amount to be run
    /// @see log_finish
    virtual void        log_start( std::ostream& os, counter_t test_cases_amount ) = 0;

    /// Invoked at the end of test module execution
    ///
    /// @param[in] os   output stream to write a messages into
    /// @see log_start
    virtual void        log_finish( std::ostream& os ) = 0;

    /// Invoked when Unit Test Framework build information is requested
    ///
    /// @param[in] os   output stream to write a messages into
    virtual void        log_build_info( std::ostream& os ) = 0;
    // @}

    // @name Test unit start/finish

    /// Invoked when test unit starts (either test suite or test case)
    ///
    /// @param[in] os   output stream to write a messages into
    /// @param[in] tu   test unit being started
    /// @see test_unit_finish
    virtual void        test_unit_start( std::ostream& os, test_unit const& tu ) = 0;

    /// Invoked when test unit finishes
    ///
    /// @param[in] os   output stream to write a messages into
    /// @param[in] tu   test unit being finished
    /// @param[in] elapsed time in microseconds spend executing this test unit
    /// @see test_unit_start
    virtual void        test_unit_finish( std::ostream& os, test_unit const& tu, unsigned long elapsed ) = 0;

    /// Invoked if test unit skipped for any reason
    ///
    /// @param[in] os   output stream to write a messages into
    /// @param[in] tu   skipped test unit
    /// @param[in] reason explanation why was it skipped
    virtual void        test_unit_skipped( std::ostream& os, test_unit const& tu, const_string /* reason */)
    {
        test_unit_skipped( os, tu );
    }

    /// Deprecated version of this interface
    virtual void        test_unit_skipped( std::ostream& /* os */, test_unit const& /* tu */) {}

    /// Invoked when a test unit is aborted
    virtual void        test_unit_aborted( std::ostream& /* os */, test_unit const& /* tu */) {}

    // @}

    // @name Uncaught exception report

    /// Invoked when Unit Test Framework detects uncaught exception
    ///
    /// The framwork calls this function when an uncaught exception it detected.
    /// This call is followed by context information:
    /// - one call to @c entry_context_start,
    /// - as many calls to @c log_entry_context as there are context entries
    /// - one call to @c entry_context_finish
    ///
    /// The logging of the exception information is finilized by a call to @c log_exception_finish.
    ///
    /// @param[in] os   output stream to write a messages into
    /// @param[in] lcd  information about the last checkpoint before the exception was triggered
    /// @param[in] ex   information about the caught exception
    /// @see log_exception_finish
    virtual void        log_exception_start( std::ostream& os, log_checkpoint_data const& lcd, execution_exception const& ex ) = 0;

    /// Invoked when Unit Test Framework detects uncaught exception
    ///
    /// Call to this function finishes uncaught exception report.
    /// @param[in] os   output stream to write a messages into
    /// @see log_exception_start
    virtual void        log_exception_finish( std::ostream& os ) = 0;
    // @}

    // @name Regular log entry

    /// Invoked by Unit Test Framework to start new log entry

    /// Call to this function starts new log entry. It is followed by series of log_entry_value calls and finally call to log_entry_finish.
    /// A log entry may consist of one or more values being reported. Some of these values will be plain strings, while others can be complicated
    /// expressions in a form of "lazy" expression template lazy_ostream.
    /// @param[in] os   output stream to write a messages into
    /// @param[in] led  log entry attributes
    /// @param[in] let  log entry type log_entry_finish
    /// @see log_entry_value, log_entry_finish
    ///
    /// @note call to this function may happen before any call to test_unit_start or all calls to test_unit_finish as the
    /// framework might log errors raised during global initialization/shutdown.
    virtual void        log_entry_start( std::ostream& os, log_entry_data const& led, log_entry_types let ) = 0;

    /// Invoked by Unit Test Framework to report a log entry content
    ///
    /// This is one of two overloaded methods to report log entry content. This one is used to report plain string value.
    /// @param[in] os   output stream to write a messages into.
    /// @param[in] value log entry string value
    /// @see log_entry_start, log_entry_finish
    virtual void        log_entry_value( std::ostream& os, const_string value ) = 0;

    /// Invoked by Unit Test Framework to report a log entry content

    /// This is one of two overloaded methods to report log entry content. This one is used to report some complicated expression passed as
    /// an expression template lazy_ostream. In most cases default implementation provided by the framework should work as is (it just converts
    /// the lazy expression into a string.
    /// @param[in] os   output stream to write a messages into
    /// @param[in] value log entry "lazy" value
    /// @see log_entry_start, log_entry_finish
    virtual void        log_entry_value( std::ostream& os, lazy_ostream const& value ); // there is a default impl

    /// Invoked by Unit Test Framework to finish a log entry report

    /// @param[in] os   output stream to write a messages into
    /// @see log_entry_start, log_entry_start
    virtual void        log_entry_finish( std::ostream& os ) = 0;
    // @}

    // @name Log entry context report

    /// Invoked by Unit Test Framework to start log entry context report

    /// Unit Test Framework logs for failed assertions and uncaught exceptions context if one was defined by a test module.
    /// Context consists of multiple "scopes" identified by description messages assigned by the test module using
    /// BOOST_TEST_INFO/BOOST_TEST_CONTEXT statements.
    /// @param[in] os   output stream to write a messages into
    /// @param[in] l    entry log_level, to be used to fine tune the message
    /// @see log_entry_context, entry_context_finish
    virtual void        entry_context_start( std::ostream& os, log_level l ) = 0;

    /// Invoked by Unit Test Framework to report log entry context "scope" description

    /// Each "scope" description is reported by separate call to log_entry_context.
    /// @param[in] os   output stream to write a messages into
    /// @param[in] value  context "scope" description
    /// @see log_entry_start, entry_context_finish
    virtual void        log_entry_context( std::ostream& os, const_string value ) = 0;

    /// Invoked by Unit Test Framework to finish log entry context report

    /// @param[in] os   output stream to write a messages into
    /// @see log_entry_start, entry_context_context
    virtual void        entry_context_finish( std::ostream& os ) = 0;
    // @}

    // @name Log level management

    /// Sets the log level of the logger/formatter
    ///
    /// Some loggers need to manage the log level by their own. This
    /// member function let the implementation decide of that.
    /// @par Since Boost 1.62
    virtual void        set_log_level(log_level new_log_level);

    /// Returns the log level of the logger/formatter
    /// @par Since Boost 1.62
    virtual log_level   get_log_level() const;
    // @}

    // @name Stream management

    /// Returns a default stream for this logger.
    ///
    /// The returned string describes the stream as if it was passed from
    /// the command line @c "--log_sink" parameter. With that regards, @b stdout and @b stderr
    /// have special meaning indicating the standard output or error stream respectively.
    ///
    /// @par Since Boost 1.62
    virtual std::string  get_default_stream_description() const
    {
        return "stdout";
    }

    // @}

protected:
    log_level           m_log_level;

};

} // namespace unit_test
} // namespace boost

//____________________________________________________________________________//

#include <boost/test/detail/enable_warnings.hpp>

#endif // BOOST_TEST_UNIT_TEST_LOG_FORMATTER_HPP_071894GER


code:
# /* **************************************************************************
#  *                                                                          *
#  *     (C) Copyright Paul Mensonides 2002.
#  *     Distributed under 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)
#  *                                                                          *
#  ************************************************************************** */
#
# /* See http://www.boost.org for most recent version. */
#
# ifndef BOOST_PREPROCESSOR_FACILITIES_INTERCEPT_HPP
# define BOOST_PREPROCESSOR_FACILITIES_INTERCEPT_HPP
#
# /* BOOST_PP_INTERCEPT */
#
# define BOOST_PP_INTERCEPT BOOST_PP_INTERCEPT_
#
# define BOOST_PP_INTERCEPT_0
# define BOOST_PP_INTERCEPT_1
# define BOOST_PP_INTERCEPT_2
# define BOOST_PP_INTERCEPT_3
# define BOOST_PP_INTERCEPT_4
# define BOOST_PP_INTERCEPT_5
# define BOOST_PP_INTERCEPT_6
# define BOOST_PP_INTERCEPT_7
# define BOOST_PP_INTERCEPT_8
# define BOOST_PP_INTERCEPT_9
# define BOOST_PP_INTERCEPT_10
# define BOOST_PP_INTERCEPT_11
# define BOOST_PP_INTERCEPT_12
# define BOOST_PP_INTERCEPT_13
# define BOOST_PP_INTERCEPT_14
# define BOOST_PP_INTERCEPT_15
# define BOOST_PP_INTERCEPT_16
# define BOOST_PP_INTERCEPT_17
# define BOOST_PP_INTERCEPT_18
# define BOOST_PP_INTERCEPT_19
# define BOOST_PP_INTERCEPT_20
# define BOOST_PP_INTERCEPT_21
# define BOOST_PP_INTERCEPT_22
# define BOOST_PP_INTERCEPT_23
# define BOOST_PP_INTERCEPT_24
# define BOOST_PP_INTERCEPT_25
# define BOOST_PP_INTERCEPT_26
# define BOOST_PP_INTERCEPT_27
# define BOOST_PP_INTERCEPT_28
# define BOOST_PP_INTERCEPT_29
# define BOOST_PP_INTERCEPT_30
# define BOOST_PP_INTERCEPT_31
# define BOOST_PP_INTERCEPT_32
# define BOOST_PP_INTERCEPT_33
# define BOOST_PP_INTERCEPT_34
# define BOOST_PP_INTERCEPT_35
# define BOOST_PP_INTERCEPT_36
# define BOOST_PP_INTERCEPT_37
# define BOOST_PP_INTERCEPT_38
# define BOOST_PP_INTERCEPT_39
# define BOOST_PP_INTERCEPT_40
# define BOOST_PP_INTERCEPT_41
# define BOOST_PP_INTERCEPT_42
# define BOOST_PP_INTERCEPT_43
# define BOOST_PP_INTERCEPT_44
# define BOOST_PP_INTERCEPT_45
# define BOOST_PP_INTERCEPT_46
# define BOOST_PP_INTERCEPT_47
# define BOOST_PP_INTERCEPT_48
# define BOOST_PP_INTERCEPT_49
# define BOOST_PP_INTERCEPT_50
# define BOOST_PP_INTERCEPT_51
# define BOOST_PP_INTERCEPT_52
# define BOOST_PP_INTERCEPT_53
# define BOOST_PP_INTERCEPT_54
# define BOOST_PP_INTERCEPT_55
# define BOOST_PP_INTERCEPT_56
# define BOOST_PP_INTERCEPT_57
# define BOOST_PP_INTERCEPT_58
# define BOOST_PP_INTERCEPT_59
# define BOOST_PP_INTERCEPT_60
# define BOOST_PP_INTERCEPT_61
# define BOOST_PP_INTERCEPT_62
# define BOOST_PP_INTERCEPT_63
# define BOOST_PP_INTERCEPT_64
# define BOOST_PP_INTERCEPT_65
# define BOOST_PP_INTERCEPT_66
# define BOOST_PP_INTERCEPT_67
# define BOOST_PP_INTERCEPT_68
# define BOOST_PP_INTERCEPT_69
# define BOOST_PP_INTERCEPT_70
# define BOOST_PP_INTERCEPT_71
# define BOOST_PP_INTERCEPT_72
# define BOOST_PP_INTERCEPT_73
# define BOOST_PP_INTERCEPT_74
# define BOOST_PP_INTERCEPT_75
# define BOOST_PP_INTERCEPT_76
# define BOOST_PP_INTERCEPT_77
# define BOOST_PP_INTERCEPT_78
# define BOOST_PP_INTERCEPT_79
# define BOOST_PP_INTERCEPT_80
# define BOOST_PP_INTERCEPT_81
# define BOOST_PP_INTERCEPT_82
# define BOOST_PP_INTERCEPT_83
# define BOOST_PP_INTERCEPT_84
# define BOOST_PP_INTERCEPT_85
# define BOOST_PP_INTERCEPT_86
# define BOOST_PP_INTERCEPT_87
# define BOOST_PP_INTERCEPT_88
# define BOOST_PP_INTERCEPT_89
# define BOOST_PP_INTERCEPT_90
# define BOOST_PP_INTERCEPT_91
# define BOOST_PP_INTERCEPT_92
# define BOOST_PP_INTERCEPT_93
# define BOOST_PP_INTERCEPT_94
# define BOOST_PP_INTERCEPT_95
# define BOOST_PP_INTERCEPT_96
# define BOOST_PP_INTERCEPT_97
# define BOOST_PP_INTERCEPT_98
# define BOOST_PP_INTERCEPT_99
# define BOOST_PP_INTERCEPT_100
# define BOOST_PP_INTERCEPT_101
# define BOOST_PP_INTERCEPT_102
# define BOOST_PP_INTERCEPT_103
# define BOOST_PP_INTERCEPT_104
# define BOOST_PP_INTERCEPT_105
# define BOOST_PP_INTERCEPT_106
# define BOOST_PP_INTERCEPT_107
# define BOOST_PP_INTERCEPT_108
# define BOOST_PP_INTERCEPT_109
# define BOOST_PP_INTERCEPT_110
# define BOOST_PP_INTERCEPT_111
# define BOOST_PP_INTERCEPT_112
# define BOOST_PP_INTERCEPT_113
# define BOOST_PP_INTERCEPT_114
# define BOOST_PP_INTERCEPT_115
# define BOOST_PP_INTERCEPT_116
# define BOOST_PP_INTERCEPT_117
# define BOOST_PP_INTERCEPT_118
# define BOOST_PP_INTERCEPT_119
# define BOOST_PP_INTERCEPT_120
# define BOOST_PP_INTERCEPT_121
# define BOOST_PP_INTERCEPT_122
# define BOOST_PP_INTERCEPT_123
# define BOOST_PP_INTERCEPT_124
# define BOOST_PP_INTERCEPT_125
# define BOOST_PP_INTERCEPT_126
# define BOOST_PP_INTERCEPT_127
# define BOOST_PP_INTERCEPT_128
# define BOOST_PP_INTERCEPT_129
# define BOOST_PP_INTERCEPT_130
# define BOOST_PP_INTERCEPT_131
# define BOOST_PP_INTERCEPT_132
# define BOOST_PP_INTERCEPT_133
# define BOOST_PP_INTERCEPT_134
# define BOOST_PP_INTERCEPT_135
# define BOOST_PP_INTERCEPT_136
# define BOOST_PP_INTERCEPT_137
# define BOOST_PP_INTERCEPT_138
# define BOOST_PP_INTERCEPT_139
# define BOOST_PP_INTERCEPT_140
# define BOOST_PP_INTERCEPT_141
# define BOOST_PP_INTERCEPT_142
# define BOOST_PP_INTERCEPT_143
# define BOOST_PP_INTERCEPT_144
# define BOOST_PP_INTERCEPT_145
# define BOOST_PP_INTERCEPT_146
# define BOOST_PP_INTERCEPT_147
# define BOOST_PP_INTERCEPT_148
# define BOOST_PP_INTERCEPT_149
# define BOOST_PP_INTERCEPT_150
# define BOOST_PP_INTERCEPT_151
# define BOOST_PP_INTERCEPT_152
# define BOOST_PP_INTERCEPT_153
# define BOOST_PP_INTERCEPT_154
# define BOOST_PP_INTERCEPT_155
# define BOOST_PP_INTERCEPT_156
# define BOOST_PP_INTERCEPT_157
# define BOOST_PP_INTERCEPT_158
# define BOOST_PP_INTERCEPT_159
# define BOOST_PP_INTERCEPT_160
# define BOOST_PP_INTERCEPT_161
# define BOOST_PP_INTERCEPT_162
# define BOOST_PP_INTERCEPT_163
# define BOOST_PP_INTERCEPT_164
# define BOOST_PP_INTERCEPT_165
# define BOOST_PP_INTERCEPT_166
# define BOOST_PP_INTERCEPT_167
# define BOOST_PP_INTERCEPT_168
# define BOOST_PP_INTERCEPT_169
# define BOOST_PP_INTERCEPT_170
# define BOOST_PP_INTERCEPT_171
# define BOOST_PP_INTERCEPT_172
# define BOOST_PP_INTERCEPT_173
# define BOOST_PP_INTERCEPT_174
# define BOOST_PP_INTERCEPT_175
# define BOOST_PP_INTERCEPT_176
# define BOOST_PP_INTERCEPT_177
# define BOOST_PP_INTERCEPT_178
# define BOOST_PP_INTERCEPT_179
# define BOOST_PP_INTERCEPT_180
# define BOOST_PP_INTERCEPT_181
# define BOOST_PP_INTERCEPT_182
# define BOOST_PP_INTERCEPT_183
# define BOOST_PP_INTERCEPT_184
# define BOOST_PP_INTERCEPT_185
# define BOOST_PP_INTERCEPT_186
# define BOOST_PP_INTERCEPT_187
# define BOOST_PP_INTERCEPT_188
# define BOOST_PP_INTERCEPT_189
# define BOOST_PP_INTERCEPT_190
# define BOOST_PP_INTERCEPT_191
# define BOOST_PP_INTERCEPT_192
# define BOOST_PP_INTERCEPT_193
# define BOOST_PP_INTERCEPT_194
# define BOOST_PP_INTERCEPT_195
# define BOOST_PP_INTERCEPT_196
# define BOOST_PP_INTERCEPT_197
# define BOOST_PP_INTERCEPT_198
# define BOOST_PP_INTERCEPT_199
# define BOOST_PP_INTERCEPT_200
# define BOOST_PP_INTERCEPT_201
# define BOOST_PP_INTERCEPT_202
# define BOOST_PP_INTERCEPT_203
# define BOOST_PP_INTERCEPT_204
# define BOOST_PP_INTERCEPT_205
# define BOOST_PP_INTERCEPT_206
# define BOOST_PP_INTERCEPT_207
# define BOOST_PP_INTERCEPT_208
# define BOOST_PP_INTERCEPT_209
# define BOOST_PP_INTERCEPT_210
# define BOOST_PP_INTERCEPT_211
# define BOOST_PP_INTERCEPT_212
# define BOOST_PP_INTERCEPT_213
# define BOOST_PP_INTERCEPT_214
# define BOOST_PP_INTERCEPT_215
# define BOOST_PP_INTERCEPT_216
# define BOOST_PP_INTERCEPT_217
# define BOOST_PP_INTERCEPT_218
# define BOOST_PP_INTERCEPT_219
# define BOOST_PP_INTERCEPT_220
# define BOOST_PP_INTERCEPT_221
# define BOOST_PP_INTERCEPT_222
# define BOOST_PP_INTERCEPT_223
# define BOOST_PP_INTERCEPT_224
# define BOOST_PP_INTERCEPT_225
# define BOOST_PP_INTERCEPT_226
# define BOOST_PP_INTERCEPT_227
# define BOOST_PP_INTERCEPT_228
# define BOOST_PP_INTERCEPT_229
# define BOOST_PP_INTERCEPT_230
# define BOOST_PP_INTERCEPT_231
# define BOOST_PP_INTERCEPT_232
# define BOOST_PP_INTERCEPT_233
# define BOOST_PP_INTERCEPT_234
# define BOOST_PP_INTERCEPT_235
# define BOOST_PP_INTERCEPT_236
# define BOOST_PP_INTERCEPT_237
# define BOOST_PP_INTERCEPT_238
# define BOOST_PP_INTERCEPT_239
# define BOOST_PP_INTERCEPT_240
# define BOOST_PP_INTERCEPT_241
# define BOOST_PP_INTERCEPT_242
# define BOOST_PP_INTERCEPT_243
# define BOOST_PP_INTERCEPT_244
# define BOOST_PP_INTERCEPT_245
# define BOOST_PP_INTERCEPT_246
# define BOOST_PP_INTERCEPT_247
# define BOOST_PP_INTERCEPT_248
# define BOOST_PP_INTERCEPT_249
# define BOOST_PP_INTERCEPT_250
# define BOOST_PP_INTERCEPT_251
# define BOOST_PP_INTERCEPT_252
# define BOOST_PP_INTERCEPT_253
# define BOOST_PP_INTERCEPT_254
# define BOOST_PP_INTERCEPT_255
# define BOOST_PP_INTERCEPT_256
#
# endif


code:
module FFTW3
  use, intrinsic :: iso_c_binding
  include 'fftw3.f03'
end module FFTW3


code:
//  abs.hpp
//
//  (C) Copyright 2011 Vicente J. Botet Escriba
//  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_RATIO_MPL_RATIONAL_C_TAG_HPP
#define BOOST_RATIO_MPL_RATIONAL_C_TAG_HPP

#ifdef BOOST_RATIO_EXTENSIONS

#include <boost/mpl/int.hpp>

namespace boost {
namespace mpl {

struct rational_c_tag : int_<10> {};

}
}

#endif // BOOST_RATIO_EXTENSIONS
#endif  // BOOST_RATIO_MPL_RATIONAL_C_TAG_HPP


code:
//////////////////////////////////////////////////////////////////////////////
//
// (C) Copyright Ion Gaztanaga 2014-2014.
// Distributed under 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)
//
// See http://www.boost.org/libs/move for documentation.
//
//////////////////////////////////////////////////////////////////////////////

//! \file

#ifndef BOOST_MOVE_DETAIL_ITERATOR_TRAITS_HPP
#define BOOST_MOVE_DETAIL_ITERATOR_TRAITS_HPP

#ifndef BOOST_CONFIG_HPP
#  include <boost/config.hpp>
#endif
#
#if defined(BOOST_HAS_PRAGMA_ONCE)
#  pragma once
#endif

#include <cstddef>
#include <boost/move/detail/type_traits.hpp>

#include <boost/move/detail/std_ns_begin.hpp>
BOOST_MOVE_STD_NS_BEG

struct input_iterator_tag;
struct forward_iterator_tag;
struct bidirectional_iterator_tag;
struct random_access_iterator_tag;
struct output_iterator_tag;

BOOST_MOVE_STD_NS_END
#include <boost/move/detail/std_ns_end.hpp>

namespace boost{  namespace movelib{

template<class Iterator>
struct iterator_traits
{
   typedef typename Iterator::difference_type   difference_type;
   typedef typename Iterator::value_type        value_type;
   typedef typename Iterator::pointer           pointer;
   typedef typename Iterator::reference         reference;
   typedef typename Iterator::iterator_category iterator_category;
   typedef typename boost::move_detail::make_unsigned<difference_type>::type size_type;
};

template<class T>
struct iterator_traits<T*>
{
   typedef std::ptrdiff_t                    difference_type;
   typedef T                                 value_type;
   typedef T*                                pointer;
   typedef T&                                reference;
   typedef std::random_access_iterator_tag   iterator_category;
   typedef typename boost::move_detail::make_unsigned<difference_type>::type size_type;
};

template<class T>
struct iterator_traits<const T*>
{
   typedef std::ptrdiff_t                    difference_type;
   typedef T                                 value_type;
   typedef const T*                          pointer;
   typedef const T&                          reference;
   typedef std::random_access_iterator_tag   iterator_category;
   typedef typename boost::move_detail::make_unsigned<difference_type>::type size_type;
};

}} //namespace boost {  namespace movelib{

#endif //#ifndef BOOST_MOVE_DETAIL_ITERATOR_TRAITS_HPP


code:
subroutine getfc2w(c,csync,npeaks,fs,fc1,fpks)

  include 'wsprlf_params.f90'

  complex c(0:NZ-1)                     !Complex waveform
  complex cs(0:NZ-1)                    !For computing spectrum
  complex csync(0:NZ-1)                 !Sync symbols only, from cbb
  real a(5)
  real freqs(413),sp2(413),fpks(npeaks)
  integer pkloc(1)

  df=fs/NZ
  baud=fs/NSPS
  a(1)=-fc1
  a(2:5)=0.
  call twkfreq1(c,NZ,fs,a,cs)         !Mix down by fc1

! Filter, square, then FFT to get refined carrier frequency fc2.
  call four2a(cs,NZ,1,-1,1)          !To freq domain

  ia=nint(0.75*baud/df)
  cs(ia:NZ-1-ia)=0.                  !Save only freqs around fc1
  call four2a(cs,NZ,1,1,1)           !Back to time domain
  cs=cs/NZ
  cs=cs*cs                           !Square the data
  call four2a(cs,NZ,1,-1,1)          !Compute squared spectrum

! Find two peaks separated by baud
  pmax=0.
  fc2=0.
  ja=nint(0.3*baud/df)
  k=1
  do j=-ja,ja
     f2=j*df
     ia=nint((f2-0.5*baud)/df)
     if(ia.lt.0) ia=ia+NZ
     ib=nint((f2+0.5*baud)/df)
     p=real(cs(ia))**2 + aimag(cs(ia))**2 +                        &
          real(cs(ib))**2 + aimag(cs(ib))**2
     if(p.gt.pmax) then
        pmax=p
        fc2=0.5*f2
     endif
     freqs(k)=0.5*f2
     sp2(k)=p
     k=k+1
!           write(52,1200) f2,p,db(p)
!1200       format(f10.3,2f15.3)
  enddo

  do i=1,npeaks
    pkloc=maxloc(sp2)
    ipk=pkloc(1)
    fpks(i)=freqs(ipk)
    ipk0=max(1,ipk-1)
    ipk1=min(413,ipk+1)
!    ipk0=ipk
!    ipk1=ipk
    sp2(ipk0:ipk1)=0.0
!write(*,*) i,fpks(i),fc2
  enddo

  a(1)=-fc1
  a(2:5)=0.
  call twkfreq1(c,NZ,fs,a,cs)         !Mix down by fc1
  cs=cs*conjg(csync)
  call four2a(cs,NZ,1,-1,1)          !To freq domain
  pmax=0.
  do i=0,NZ-1
     f=i*df
     if(i.gt.NZ/2) f=(i-NZ)*df
     p=real(cs(i))**2 + aimag(cs(i))**2
!     write(51,3001) f,p,db(p)
!3001 format(f10.3,e12.3,f10.3)
     if(p.gt.pmax) then
        pmax=p
        fc3=f
     endif
  enddo

  return
end subroutine getfc2w