code:
///////////////////////////////////////////////////////////////
//  Copyright 2012 John Maddock. 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_

#ifndef BOOST_MP_DIGITS_HPP
#define BOOST_MP_DIGITS_HPP

namespace boost{ namespace multiprecision{ namespace detail{

inline unsigned long digits10_2_2(unsigned long d10)
{
   return (d10 * 1000uL) / 301uL + ((d10 * 1000uL) % 301 ? 2u : 1u);
}

inline unsigned long digits2_2_10(unsigned long d2)
{
   return (d2 * 301uL) / 1000uL;
}

}}} // namespaces

#endif

license: mit


code:
//  (C) Copyright John Maddock 2003.
//  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_CONFIG_REQUIRES_THREADS_HPP
#define BOOST_CONFIG_REQUIRES_THREADS_HPP

#ifndef BOOST_CONFIG_HPP
#  include <boost/config.hpp>
#endif

#if defined(BOOST_DISABLE_THREADS)

//
// special case to handle versions of gcc which don't currently support threads:
//
#if defined(__GNUC__) && ((__GNUC__ < 3) || (__GNUC_MINOR__ <= 3) || !defined(BOOST_STRICT_CONFIG))
//
// this is checked up to gcc 3.3:
//
#if defined(__sgi) || defined(__hpux)
#  error "Multi-threaded programs are not supported by gcc on HPUX or Irix (last checked with gcc 3.3)"
#endif

#endif

#  error "Threading support unavaliable: it has been explicitly disabled with BOOST_DISABLE_THREADS"

#elif !defined(BOOST_HAS_THREADS)

# if defined __COMO__
//  Comeau C++
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -D_MT (Windows) or -D_REENTRANT (Unix)"

#elif defined(__INTEL_COMPILER) || defined(__ICL) || defined(__ICC) || defined(__ECC)
//  Intel
#ifdef _WIN32
#  error "Compiler threading support is not turned on. Please set the correct command line options for threading: either /MT /MTd /MD or /MDd"
#else
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -openmp"
#endif

# elif defined __GNUC__
//  GNU C++:
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -pthread (Linux), -pthreads (Solaris) or -mthreads (Mingw32)"

#elif defined __sgi
//  SGI MIPSpro C++
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -D_SGI_MP_SOURCE"

#elif defined __DECCXX
//  Compaq Tru64 Unix cxx
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -pthread"

#elif defined __BORLANDC__
//  Borland
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -tWM"

#elif defined  __MWERKS__
//  Metrowerks CodeWarrior
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: either -runtime sm, -runtime smd, -runtime dm, or -runtime dmd"

#elif defined  __SUNPRO_CC
//  Sun Workshop Compiler C++
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -mt"

#elif defined __HP_aCC
//  HP aCC
#   error "Compiler threading support is not turned on. Please set the correct command line options for threading: -mt"

#elif defined(__IBMCPP__)
//  IBM Visual Age
#   error "Compiler threading support is not turned on. Please compile the code with the xlC_r compiler"

#elif defined _MSC_VER
//  Microsoft Visual C++
//
//  Must remain the last #elif since some other vendors (Metrowerks, for
//  example) also #define _MSC_VER
#  error "Compiler threading support is not turned on. Please set the correct command line options for threading: either /MT /MTd /MD or /MDd"

#else

#  error "Compiler threading support is not turned on.  Please consult your compiler's documentation for the appropriate options to use"

#endif // compilers

#endif // BOOST_HAS_THREADS

#endif // BOOST_CONFIG_REQUIRES_THREADS_HPP


code:
//  Copyright John Maddock 2006-7, 2013-14.
//  Copyright Paul A. Bristow 2007, 2013-14.
//  Copyright Nikhar Agrawal 2013-14
//  Copyright Christopher Kormanyos 2013-14

//  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_SF_GAMMA_HPP
#define BOOST_MATH_SF_GAMMA_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/config.hpp>
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/fraction.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/special_functions/powm1.hpp>
#include <boost/math/special_functions/sqrt1pm1.hpp>
#include <boost/math/special_functions/lanczos.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/detail/igamma_large.hpp>
#include <boost/math/special_functions/detail/unchecked_factorial.hpp>
#include <boost/math/special_functions/detail/lgamma_small.hpp>
#include <boost/math/special_functions/bernoulli.hpp>
#include <boost/math/special_functions/zeta.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/assert.hpp>
#include <boost/mpl/greater.hpp>
#include <boost/mpl/equal_to.hpp>
#include <boost/mpl/greater.hpp>

#include <boost/config/no_tr1/cmath.hpp>
#include <algorithm>

#ifdef BOOST_MSVC
# pragma warning(push)
# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
# pragma warning(disable: 4127) // conditional expression is constant.
# pragma warning(disable: 4100) // unreferenced formal parameter.
// Several variables made comments,
// but some difficulty as whether referenced on not may depend on macro values.
// So to be safe, 4100 warnings suppressed.
// TODO - revisit this?
#endif

namespace boost{ namespace math{

namespace detail{

template <class T>
inline bool is_odd(T v, const boost::true_type&)
{
   int i = static_cast<int>(v);
   return i&1;
}
template <class T>
inline bool is_odd(T v, const boost::false_type&)
{
   // Oh dear can't cast T to int!
   BOOST_MATH_STD_USING
   T modulus = v - 2 * floor(v/2);
   return static_cast<bool>(modulus != 0);
}
template <class T>
inline bool is_odd(T v)
{
   return is_odd(v, ::boost::is_convertible<T, int>());
}

template <class T>
T sinpx(T z)
{
   // Ad hoc function calculates x * sin(pi * x),
   // taking extra care near when x is near a whole number.
   BOOST_MATH_STD_USING
   int sign = 1;
   if(z < 0)
   {
      z = -z;
   }
   T fl = floor(z);
   T dist;
   if(is_odd(fl))
   {
      fl += 1;
      dist = fl - z;
      sign = -sign;
   }
   else
   {
      dist = z - fl;
   }
   BOOST_ASSERT(fl >= 0);
   if(dist > 0.5)
      dist = 1 - dist;
   T result = sin(dist*boost::math::constants::pi<T>());
   return sign*z*result;
} // template <class T> T sinpx(T z)
//
// tgamma(z), with Lanczos support:
//
template <class T, class Policy, class Lanczos>
T gamma_imp(T z, const Policy& pol, const Lanczos& l)
{
   BOOST_MATH_STD_USING

   T result = 1;

#ifdef BOOST_MATH_INSTRUMENT
   static bool b = false;
   if(!b)
   {
      std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
      b = true;
   }
#endif
   static const char* function = "boost::math::tgamma<%1%>(%1%)";

   if(z <= 0)
   {
      if(floor(z) == z)
         return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
      if(z <= -20)
      {
         result = gamma_imp(T(-z), pol, l) * sinpx(z);
         BOOST_MATH_INSTRUMENT_VARIABLE(result);
         if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
            return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
         result = -boost::math::constants::pi<T>() / result;
         if(result == 0)
            return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
         if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
            return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
         BOOST_MATH_INSTRUMENT_VARIABLE(result);
         return result;
      }

      // shift z to > 1:
      while(z < 0)
      {
         result /= z;
         z += 1;
      }
   }
   BOOST_MATH_INSTRUMENT_VARIABLE(result);
   if((floor(z) == z) && (z < max_factorial<T>::value))
   {
      result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
      BOOST_MATH_INSTRUMENT_VARIABLE(result);
   }
   else if (z < tools::root_epsilon<T>())
   {
      if (z < 1 / tools::max_value<T>())
         result = policies::raise_overflow_error<T>(function, 0, pol);
      result *= 1 / z - constants::euler<T>();
   }
   else
   {
      result *= Lanczos::lanczos_sum(z);
      T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
      T lzgh = log(zgh);
      BOOST_MATH_INSTRUMENT_VARIABLE(result);
      BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
      if(z * lzgh > tools::log_max_value<T>())
      {
         // we're going to overflow unless this is done with care:
         BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
         if(lzgh * z / 2 > tools::log_max_value<T>())
            return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
         T hp = pow(zgh, (z / 2) - T(0.25));
         BOOST_MATH_INSTRUMENT_VARIABLE(hp);
         result *= hp / exp(zgh);
         BOOST_MATH_INSTRUMENT_VARIABLE(result);
         if(tools::max_value<T>() / hp < result)
            return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
         result *= hp;
         BOOST_MATH_INSTRUMENT_VARIABLE(result);
      }
      else
      {
         BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
         BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
         BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
         result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
         BOOST_MATH_INSTRUMENT_VARIABLE(result);
      }
   }
   return result;
}
//
// lgamma(z) with Lanczos support:
//
template <class T, class Policy, class Lanczos>
T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
{
#ifdef BOOST_MATH_INSTRUMENT
   static bool b = false;
   if(!b)
   {
      std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
      b = true;
   }
#endif

   BOOST_MATH_STD_USING

   static const char* function = "boost::math::lgamma<%1%>(%1%)";

   T result = 0;
   int sresult = 1;
   if(z <= -tools::root_epsilon<T>())
   {
      // reflection formula:
      if(floor(z) == z)
         return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);

      T t = sinpx(z);
      z = -z;
      if(t < 0)
      {
         t = -t;
      }
      else
      {
         sresult = -sresult;
      }
      result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
   }
   else if (z < tools::root_epsilon<T>())
   {
      if (0 == z)
         return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
      if (fabs(z) < 1 / tools::max_value<T>())
         result = -log(fabs(z));
      else
         result = log(fabs(1 / z - constants::euler<T>()));
      if (z < 0)
         sresult = -1;
   }
   else if(z < 15)
   {
      typedef typename policies::precision<T, Policy>::type precision_type;
      typedef typename mpl::if_<
         mpl::and_<
            mpl::less_equal<precision_type, mpl::int_<64> >,
            mpl::greater<precision_type, mpl::int_<0> >
         >,
         mpl::int_<64>,
         typename mpl::if_<
            mpl::and_<
               mpl::less_equal<precision_type, mpl::int_<113> >,
               mpl::greater<precision_type, mpl::int_<0> >
            >,
            mpl::int_<113>, mpl::int_<0> >::type
          >::type tag_type;
      result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
   }
   else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
   {
      // taking the log of tgamma reduces the error, no danger of overflow here:
      result = log(gamma_imp(z, pol, l));
   }
   else
   {
      // regular evaluation:
      T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
      result = log(zgh) - 1;
      result *= z - 0.5f;
      result += log(Lanczos::lanczos_sum_expG_scaled(z));
   }

   if(sign)
      *sign = sresult;
   return result;
}

//
// Incomplete gamma functions follow:
//
template <class T>
struct upper_incomplete_gamma_fract
{
private:
   T z, a;
   int k;
public:
   typedef std::pair<T,T> result_type;

   upper_incomplete_gamma_fract(T a1, T z1)
      : z(z1-a1+1), a(a1), k(0)
   {
   }

   result_type operator()()
   {
      ++k;
      z += 2;
      return result_type(k * (a - k), z);
   }
};

template <class T>
inline T upper_gamma_fraction(T a, T z, T eps)
{
   // Multiply result by z^a * e^-z to get the full
   // upper incomplete integral.  Divide by tgamma(z)
   // to normalise.
   upper_incomplete_gamma_fract<T> f(a, z);
   return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
}

template <class T>
struct lower_incomplete_gamma_series
{
private:
   T a, z, result;
public:
   typedef T result_type;
   lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}

   T operator()()
   {
      T r = result;
      a += 1;
      result *= z/a;
      return r;
   }
};

template <class T, class Policy>
inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
{
   // Multiply result by ((z^a) * (e^-z) / a) to get the full
   // lower incomplete integral. Then divide by tgamma(a)
   // to get the normalised value.
   lower_incomplete_gamma_series<T> s(a, z);
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
   T factor = policies::get_epsilon<T, Policy>();
   T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
   policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
   return result;
}

//
// Fully generic tgamma and lgamma use Stirling's approximation
// with Bernoulli numbers.
//
template<class T>
std::size_t highest_bernoulli_index()
{
   const float digits10_of_type = (std::numeric_limits<T>::is_specialized
                                      ? static_cast<float>(std::numeric_limits<T>::digits10)
                                      : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));

   // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
   return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
}

template<class T>
T minimum_argument_for_bernoulli_recursion()
{
   const float digits10_of_type = (std::numeric_limits<T>::is_specialized
                                      ? static_cast<float>(std::numeric_limits<T>::digits10)
                                      : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));

   return T(digits10_of_type * 1.7F);
}

// Forward declaration of the lgamma_imp template specialization.
template <class T, class Policy>
T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);

template <class T, class Policy>
T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
{
   BOOST_MATH_STD_USING

   static const char* function = "boost::math::tgamma<%1%>(%1%)";

   // Check if the argument of tgamma is identically zero.
   const bool is_at_zero = (z == 0);

   if(is_at_zero)
      return policies::raise_domain_error<T>(function, "Evaluation of tgamma at zero %1%.", z, pol);

   const bool b_neg = (z < 0);

   const bool floor_of_z_is_equal_to_z = (floor(z) == z);

   // Special case handling of small factorials:
   if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
   {
      return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
   }

   // Make a local, unsigned copy of the input argument.
   T zz((!b_neg) ? z : -z);

   // Special case for ultra-small z:
   if(zz < tools::cbrt_epsilon<T>())
   {
      const T a0(1);
      const T a1(boost::math::constants::euler<T>());
      const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
      const T a2((six_euler_squared -  boost::math::constants::pi_sqr<T>()) / 12);

      const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);

      return 1 / inverse_tgamma_series;
   }

   // Scale the argument up for the calculation of lgamma,
   // and use downward recursion later for the final result.
   const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();

   int n_recur;

   if(zz < min_arg_for_recursion)
   {
      n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;

      zz += n_recur;
   }
   else
   {
      n_recur = 0;
   }

   const T log_gamma_value = lgamma_imp(zz, pol, lanczos::undefined_lanczos());

   if(log_gamma_value > tools::log_max_value<T>())
      return policies::raise_overflow_error<T>(function, 0, pol);

   T gamma_value = exp(log_gamma_value);

   // Rescale the result using downward recursion if necessary.
   if(n_recur)
   {
      // The order of divides is important, if we keep subtracting 1 from zz
      // we DO NOT get back to z (cancellation error).  Further if z < epsilon
      // we would end up dividing by zero.  Also in order to prevent spurious
      // overflow with the first division, we must save dividing by |z| till last,
      // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
      zz = fabs(z) + 1;
      for(int k = 1; k < n_recur; ++k)
      {
         gamma_value /= zz;
         zz += 1;
      }
      gamma_value /= fabs(z);
   }

   // Return the result, accounting for possible negative arguments.
   if(b_neg)
   {
      // Provide special error analysis for:
      // * arguments in the neighborhood of a negative integer
      // * arguments exactly equal to a negative integer.

      // Check if the argument of tgamma is exactly equal to a negative integer.
      if(floor_of_z_is_equal_to_z)
         return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);

      gamma_value *= sinpx(z);

      BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);

      const bool result_is_too_large_to_represent = (   (abs(gamma_value) < 1)
                                                     && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));

      if(result_is_too_large_to_represent)
         return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);

      gamma_value = -boost::math::constants::pi<T>() / gamma_value;
      BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);

      if(gamma_value == 0)
         return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);

      if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
         return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
   }

   return gamma_value;
}

template <class T, class Policy>
inline T log_gamma_near_1(const T& z, Policy const& pol)
{
   //
   // This is for the multiprecision case where there is
   // no lanczos support...
   //
   BOOST_MATH_STD_USING // ADL of std names

   BOOST_ASSERT(fabs(z) < 1);

   T result = -constants::euler<T>() * z;

   T power_term = z * z;
   T term;
   unsigned j = 0;

   do
   {
      term = boost::math::zeta<T>(j + 2, pol) * power_term / (j + 2);
      if(j & 1)
         result -= term;
      else
         result += term;
      power_term *= z;
      ++j;
   } while(fabs(result) * tools::epsilon<T>() < fabs(term));

   return result;
}

template <class T, class Policy>
T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
{
   BOOST_MATH_STD_USING

   static const char* function = "boost::math::lgamma<%1%>(%1%)";

   // Check if the argument of lgamma is identically zero.
   const bool is_at_zero = (z == 0);

   if(is_at_zero)
      return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);

   const bool b_neg = (z < 0);

   const bool floor_of_z_is_equal_to_z = (floor(z) == z);

   // Special case handling of small factorials:
   if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
   {
      return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
   }

   // Make a local, unsigned copy of the input argument.
   T zz((!b_neg) ? z : -z);

   const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();

   T log_gamma_value;

   if (zz < min_arg_for_recursion)
   {
      // Here we simply take the logarithm of tgamma(). This is somewhat
      // inefficient, but simple. The rationale is that the argument here
      // is relatively small and overflow is not expected to be likely.
      if(fabs(z - 1) < 0.25)
      {
         return log_gamma_near_1(T(zz - 1), pol);
      }
      else if(fabs(z - 2) < 0.25)
      {
         return log_gamma_near_1(T(zz - 2), pol) + log(zz - 1);
      }
      else if (z > -tools::root_epsilon<T>())
      {
         // Reflection formula may fail if z is very close to zero, let the series
         // expansion for tgamma close to zero do the work:
         log_gamma_value = log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
         if (sign)
         {
             *sign = z < 0 ? -1 : 1;
         }
         return log_gamma_value;
      }
      else
      {
         // No issue with spurious overflow in reflection formula,
         // just fall through to regular code:
         log_gamma_value = log(abs(gamma_imp(zz, pol, lanczos::undefined_lanczos())));
      }
   }
   else
   {
      // Perform the Bernoulli series expansion of Stirling's approximation.

      const std::size_t number_of_bernoullis_b2n = highest_bernoulli_index<T>();

            T one_over_x_pow_two_n_minus_one = 1 / zz;
      const T one_over_x2                    = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
            T sum                            = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
      const T target_epsilon_to_break_loop   = (sum * boost::math::tools::epsilon<T>()) * T(1.0E-10F);

      for(std::size_t n = 2U; n < number_of_bernoullis_b2n; ++n)
      {
         one_over_x_pow_two_n_minus_one *= one_over_x2;

         const std::size_t n2 = static_cast<std::size_t>(n * 2U);

         const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));

         if((n >= 8U) && (abs(term) < target_epsilon_to_break_loop))
         {
            // We have reached the desired precision in Stirling's expansion.
            // Adding additional terms to the sum of this divergent asymptotic
            // expansion will not improve the result.

            // Break from the loop.
            break;
         }

         sum += term;
      }

      // Complete Stirling's approximation.
      const T half_ln_two_pi = log(boost::math::constants::two_pi<T>()) / 2;

      log_gamma_value = ((((zz - boost::math::constants::half<T>()) * log(zz)) - zz) + half_ln_two_pi) + sum;
   }

   int sign_of_result = 1;

   if(b_neg)
   {
      // Provide special error analysis if the argument is exactly
      // equal to a negative integer.

      // Check if the argument of lgamma is exactly equal to a negative integer.
      if(floor_of_z_is_equal_to_z)
         return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);

      T t = sinpx(z);

      if(t < 0)
      {
         t = -t;
      }
      else
      {
         sign_of_result = -sign_of_result;
      }

      log_gamma_value = - log_gamma_value
                        + log(boost::math::constants::pi<T>())
                        - log(t);
   }

   if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }

   return log_gamma_value;
}

//
// This helper calculates tgamma(dz+1)-1 without cancellation errors,
// used by the upper incomplete gamma with z < 1:
//
template <class T, class Policy, class Lanczos>
T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
{
   BOOST_MATH_STD_USING

   typedef typename policies::precision<T,Policy>::type precision_type;

   typedef typename mpl::if_<
      mpl::or_<
         mpl::less_equal<precision_type, mpl::int_<0> >,
         mpl::greater<precision_type, mpl::int_<113> >
      >,
      typename mpl::if_<
         mpl::and_<is_same<Lanczos, lanczos::lanczos24m113>, mpl::greater<precision_type, mpl::int_<0> > >,
         mpl::int_<113>,
         mpl::int_<0>
      >::type,
      typename mpl::if_<
         mpl::less_equal<precision_type, mpl::int_<64> >,
         mpl::int_<64>, mpl::int_<113> >::type
       >::type tag_type;

   T result;
   if(dz < 0)
   {
      if(dz < -0.5)
      {
         // Best method is simply to subtract 1 from tgamma:
         result = boost::math::tgamma(1+dz, pol) - 1;
         BOOST_MATH_INSTRUMENT_CODE(result);
      }
      else
      {
         // Use expm1 on lgamma:
         result = boost::math::expm1(-boost::math::log1p(dz, pol)
            + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
         BOOST_MATH_INSTRUMENT_CODE(result);
      }
   }
   else
   {
      if(dz < 2)
      {
         // Use expm1 on lgamma:
         result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
         BOOST_MATH_INSTRUMENT_CODE(result);
      }
      else
      {
         // Best method is simply to subtract 1 from tgamma:
         result = boost::math::tgamma(1+dz, pol) - 1;
         BOOST_MATH_INSTRUMENT_CODE(result);
      }
   }

   return result;
}

template <class T, class Policy>
inline T tgammap1m1_imp(T z, Policy const& pol,
                 const ::boost::math::lanczos::undefined_lanczos&)
{
   BOOST_MATH_STD_USING // ADL of std names

   if(fabs(z) < 0.55)
   {
      return boost::math::expm1(log_gamma_near_1(z, pol));
   }
   return boost::math::expm1(boost::math::lgamma(1 + z, pol));
}

//
// Series representation for upper fraction when z is small:
//
template <class T>
struct small_gamma2_series
{
   typedef T result_type;

   small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}

   T operator()()
   {
      T r = result / (apn);
      result *= x;
      result /= ++n;
      apn += 1;
      return r;
   }

private:
   T result, x, apn;
   int n;
};
//
// calculate power term prefix (z^a)(e^-z) used in the non-normalised
// incomplete gammas:
//
template <class T, class Policy>
T full_igamma_prefix(T a, T z, const Policy& pol)
{
   BOOST_MATH_STD_USING

   T prefix;
   T alz = a * log(z);

   if(z >= 1)
   {
      if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
      {
         prefix = pow(z, a) * exp(-z);
      }
      else if(a >= 1)
      {
         prefix = pow(z / exp(z/a), a);
      }
      else
      {
         prefix = exp(alz - z);
      }
   }
   else
   {
      if(alz > tools::log_min_value<T>())
      {
         prefix = pow(z, a) * exp(-z);
      }
      else if(z/a < tools::log_max_value<T>())
      {
         prefix = pow(z / exp(z/a), a);
      }
      else
      {
         prefix = exp(alz - z);
      }
   }
   //
   // This error handling isn't very good: it happens after the fact
   // rather than before it...
   //
   if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
      return policies::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol);

   return prefix;
}
//
// Compute (z^a)(e^-z)/tgamma(a)
// most if the error occurs in this function:
//
template <class T, class Policy, class Lanczos>
T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
{
   BOOST_MATH_STD_USING
   T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
   T prefix;
   T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;

   if(a < 1)
   {
      //
      // We have to treat a < 1 as a special case because our Lanczos
      // approximations are optimised against the factorials with a > 1,
      // and for high precision types especially (128-bit reals for example)
      // very small values of a can give rather eroneous results for gamma
      // unless we do this:
      //
      // TODO: is this still required?  Lanczos approx should be better now?
      //
      if(z <= tools::log_min_value<T>())
      {
         // Oh dear, have to use logs, should be free of cancellation errors though:
         return exp(a * log(z) - z - lgamma_imp(a, pol, l));
      }
      else
      {
         // direct calculation, no danger of overflow as gamma(a) < 1/a
         // for small a.
         return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
      }
   }
   else if((fabs(d*d*a) <= 100) && (a > 150))
   {
      // special case for large a and a ~ z.
      prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
      prefix = exp(prefix);
   }
   else
   {
      //
      // general case.
      // direct computation is most accurate, but use various fallbacks
      // for different parts of the problem domain:
      //
      T alz = a * log(z / agh);
      T amz = a - z;
      if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
      {
         T amza = amz / a;
         if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
         {
            // compute square root of the result and then square it:
            T sq = pow(z / agh, a / 2) * exp(amz / 2);
            prefix = sq * sq;
         }
         else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
         {
            // compute the 4th root of the result then square it twice:
            T sq = pow(z / agh, a / 4) * exp(amz / 4);
            prefix = sq * sq;
            prefix *= prefix;
         }
         else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
         {
            prefix = pow((z * exp(amza)) / agh, a);
         }
         else
         {
            prefix = exp(alz + amz);
         }
      }
      else
      {
         prefix = pow(z / agh, a) * exp(amz);
      }
   }
   prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
   return prefix;
}
//
// And again, without Lanczos support:
//
template <class T, class Policy>
T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
{
   BOOST_MATH_STD_USING

   T limit = (std::max)(T(10), a);
   T sum = detail::lower_gamma_series(a, limit, pol) / a;
   sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());

   if(a < 10)
   {
      // special case for small a:
      T prefix = pow(z / 10, a);
      prefix *= exp(10-z);
      if(0 == prefix)
      {
         prefix = pow((z * exp((10-z)/a)) / 10, a);
      }
      prefix /= sum;
      return prefix;
   }

   T zoa = z / a;
   T amz = a - z;
   T alzoa = a * log(zoa);
   T prefix;
   if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
   {
      T amza = amz / a;
      if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
      {
         prefix = exp(alzoa + amz);
      }
      else
      {
         prefix = pow(zoa * exp(amza), a);
      }
   }
   else
   {
      prefix = pow(zoa, a) * exp(amz);
   }
   prefix /= sum;
   return prefix;
}
//
// Upper gamma fraction for very small a:
//
template <class T, class Policy>
inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
{
   BOOST_MATH_STD_USING  // ADL of std functions.
   //
   // Compute the full upper fraction (Q) when a is very small:
   //
   T result;
   result = boost::math::tgamma1pm1(a, pol);
   if(pgam)
      *pgam = (result + 1) / a;
   T p = boost::math::powm1(x, a, pol);
   result -= p;
   result /= a;
   detail::small_gamma2_series<T> s(a, x);
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
   p += 1;
   if(pderivative)
      *pderivative = p / (*pgam * exp(x));
   T init_value = invert ? *pgam : 0;
   result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
   policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
   if(invert)
      result = -result;
   return result;
}
//
// Upper gamma fraction for integer a:
//
template <class T, class Policy>
inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
{
   //
   // Calculates normalised Q when a is an integer:
   //
   BOOST_MATH_STD_USING
   T e = exp(-x);
   T sum = e;
   if(sum != 0)
   {
      T term = sum;
      for(unsigned n = 1; n < a; ++n)
      {
         term /= n;
         term *= x;
         sum += term;
      }
   }
   if(pderivative)
   {
      *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
   }
   return sum;
}
//
// Upper gamma fraction for half integer a:
//
template <class T, class Policy>
T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
{
   //
   // Calculates normalised Q when a is a half-integer:
   //
   BOOST_MATH_STD_USING
   T e = boost::math::erfc(sqrt(x), pol);
   if((e != 0) && (a > 1))
   {
      T term = exp(-x) / sqrt(constants::pi<T>() * x);
      term *= x;
      static const T half = T(1) / 2;
      term /= half;
      T sum = term;
      for(unsigned n = 2; n < a; ++n)
      {
         term /= n - half;
         term *= x;
         sum += term;
      }
      e += sum;
      if(p_derivative)
      {
         *p_derivative = 0;
      }
   }
   else if(p_derivative)
   {
      // We'll be dividing by x later, so calculate derivative * x:
      *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
   }
   return e;
}
//
// Main incomplete gamma entry point, handles all four incomplete gamma's:
//
template <class T, class Policy>
T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
                       const Policy& pol, T* p_derivative)
{
   static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
   if(a <= 0)
      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
   if(x < 0)
      return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);

   BOOST_MATH_STD_USING

   typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;

   T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used

   if(a >= max_factorial<T>::value && !normalised)
   {
      //
      // When we're computing the non-normalized incomplete gamma
      // and a is large the result is rather hard to compute unless
      // we use logs.  There are really two options - if x is a long
      // way from a in value then we can reliably use methods 2 and 4
      // below in logarithmic form and go straight to the result.
      // Otherwise we let the regularized gamma take the strain
      // (the result is unlikely to unerflow in the central region anyway)
      // and combine with lgamma in the hopes that we get a finite result.
      //
      if(invert && (a * 4 < x))
      {
         // This is method 4 below, done in logs:
         result = a * log(x) - x;
         if(p_derivative)
            *p_derivative = exp(result);
         result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
      }
      else if(!invert && (a > 4 * x))
      {
         // This is method 2 below, done in logs:
         result = a * log(x) - x;
         if(p_derivative)
            *p_derivative = exp(result);
         T init_value = 0;
         result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
      }
      else
      {
         result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
         if(result == 0)
         {
            if(invert)
            {
               // Try http://functions.wolfram.com/06.06.06.0039.01
               result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
               result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
               if(p_derivative)
                  *p_derivative = exp(a * log(x) - x);
            }
            else
            {
               // This is method 2 below, done in logs, we're really outside the
               // range of this method, but since the result is almost certainly
               // infinite, we should probably be OK:
               result = a * log(x) - x;
               if(p_derivative)
                  *p_derivative = exp(result);
               T init_value = 0;
               result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
            }
         }
         else
         {
            result = log(result) + boost::math::lgamma(a, pol);
         }
      }
      if(result > tools::log_max_value<T>())
         return policies::raise_overflow_error<T>(function, 0, pol);
      return exp(result);
   }

   BOOST_ASSERT((p_derivative == 0) || (normalised == true));

   bool is_int, is_half_int;
   bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
   if(is_small_a)
   {
      T fa = floor(a);
      is_int = (fa == a);
      is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
   }
   else
   {
      is_int = is_half_int = false;
   }

   int eval_method;

   if(is_int && (x > 0.6))
   {
      // calculate Q via finite sum:
      invert = !invert;
      eval_method = 0;
   }
   else if(is_half_int && (x > 0.2))
   {
      // calculate Q via finite sum for half integer a:
      invert = !invert;
      eval_method = 1;
   }
   else if((x < tools::root_epsilon<T>()) && (a > 1))
   {
      eval_method = 6;
   }
   else if(x < 0.5)
   {
      //
      // Changeover criterion chosen to give a changeover at Q ~ 0.33
      //
      if(-0.4 / log(x) < a)
      {
         eval_method = 2;
      }
      else
      {
         eval_method = 3;
      }
   }
   else if(x < 1.1)
   {
      //
      // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
      //
      if(x * 0.75f < a)
      {
         eval_method = 2;
      }
      else
      {
         eval_method = 3;
      }
   }
   else
   {
      //
      // Begin by testing whether we're in the "bad" zone
      // where the result will be near 0.5 and the usual
      // series and continued fractions are slow to converge:
      //
      bool use_temme = false;
      if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
      {
         T sigma = fabs((x-a)/a);
         if((a > 200) && (policies::digits<T, Policy>() <= 113))
         {
            //
            // This limit is chosen so that we use Temme's expansion
            // only if the result would be larger than about 10^-6.
            // Below that the regular series and continued fractions
            // converge OK, and if we use Temme's method we get increasing
            // errors from the dominant erfc term as it's (inexact) argument
            // increases in magnitude.
            //
            if(20 / a > sigma * sigma)
               use_temme = true;
         }
         else if(policies::digits<T, Policy>() <= 64)
         {
            // Note in this zone we can't use Temme's expansion for
            // types longer than an 80-bit real:
            // it would require too many terms in the polynomials.
            if(sigma < 0.4)
               use_temme = true;
         }
      }
      if(use_temme)
      {
         eval_method = 5;
      }
      else
      {
         //
         // Regular case where the result will not be too close to 0.5.
         //
         // Changeover here occurs at P ~ Q ~ 0.5
         // Note that series computation of P is about x2 faster than continued fraction
         // calculation of Q, so try and use the CF only when really necessary, especially
         // for small x.
         //
         if(x - (1 / (3 * x)) < a)
         {
            eval_method = 2;
         }
         else
         {
            eval_method = 4;
            invert = !invert;
         }
      }
   }

   switch(eval_method)
   {
   case 0:
      {
         result = finite_gamma_q(a, x, pol, p_derivative);
         if(normalised == false)
            result *= boost::math::tgamma(a, pol);
         break;
      }
   case 1:
      {
         result = finite_half_gamma_q(a, x, p_derivative, pol);
         if(normalised == false)
            result *= boost::math::tgamma(a, pol);
         if(p_derivative && (*p_derivative == 0))
            *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
         break;
      }
   case 2:
      {
         // Compute P:
         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
         if(p_derivative)
            *p_derivative = result;
         if(result != 0)
         {
            //
            // If we're going to be inverting the result then we can
            // reduce the number of series evaluations by quite
            // a few iterations if we set an initial value for the
            // series sum based on what we'll end up subtracting it from
            // at the end.
            // Have to be careful though that this optimization doesn't
            // lead to spurious numberic overflow.  Note that the
            // scary/expensive overflow checks below are more often
            // than not bypassed in practice for "sensible" input
            // values:
            //
            T init_value = 0;
            bool optimised_invert = false;
            if(invert)
            {
               init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
               if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
               {
                  init_value /= result;
                  if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
                  {
                     init_value *= -a;
                     optimised_invert = true;
                  }
                  else
                     init_value = 0;
               }
               else
                  init_value = 0;
            }
            result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
            if(optimised_invert)
            {
               invert = false;
               result = -result;
            }
         }
         break;
      }
   case 3:
      {
         // Compute Q:
         invert = !invert;
         T g;
         result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
         invert = false;
         if(normalised)
            result /= g;
         break;
      }
   case 4:
      {
         // Compute Q:
         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
         if(p_derivative)
            *p_derivative = result;
         if(result != 0)
            result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
         break;
      }
   case 5:
      {
         //
         // Use compile time dispatch to the appropriate
         // Temme asymptotic expansion.  This may be dead code
         // if T does not have numeric limits support, or has
         // too many digits for the most precise version of
         // these expansions, in that case we'll be calling
         // an empty function.
         //
         typedef typename policies::precision<T, Policy>::type precision_type;

         typedef typename mpl::if_<
            mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
            mpl::greater<precision_type, mpl::int_<113> > >,
            mpl::int_<0>,
            typename mpl::if_<
               mpl::less_equal<precision_type, mpl::int_<53> >,
               mpl::int_<53>,
               typename mpl::if_<
                  mpl::less_equal<precision_type, mpl::int_<64> >,
                  mpl::int_<64>,
                  mpl::int_<113>
               >::type
            >::type
         >::type tag_type;

         result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
         if(x >= a)
            invert = !invert;
         if(p_derivative)
            *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
         break;
      }
   case 6:
      {
         // x is so small that P is necessarily very small too,
         // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
         result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol);
         result *= 1 - a * x / (a + 1);
      }
   }

   if(normalised && (result > 1))
      result = 1;
   if(invert)
   {
      T gam = normalised ? 1 : boost::math::tgamma(a, pol);
      result = gam - result;
   }
   if(p_derivative)
   {
      //
      // Need to convert prefix term to derivative:
      //
      if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
      {
         // overflow, just return an arbitrarily large value:
         *p_derivative = tools::max_value<T>() / 2;
      }

      *p_derivative /= x;
   }

   return result;
}

//
// Ratios of two gamma functions:
//
template <class T, class Policy, class Lanczos>
T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
{
   BOOST_MATH_STD_USING
   if(z < tools::epsilon<T>())
   {
      //
      // We get spurious numeric overflow unless we're very careful, this
      // can occur either inside Lanczos::lanczos_sum(z) or in the
      // final combination of terms, to avoid this, split the product up
      // into 2 (or 3) parts:
      //
      // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
      //    z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
      //
      if(boost::math::max_factorial<T>::value < delta)
      {
         T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
         ratio *= z;
         ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
         return 1 / ratio;
      }
      else
      {
         return 1 / (z * boost::math::tgamma(z + delta, pol));
      }
   }
   T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>());
   T result;
   if(z + delta == z)
   {
      if(fabs(delta) < 10)
         result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
      else
         result = 1;
   }
   else
   {
      if(fabs(delta) < 10)
      {
         result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
      }
      else
      {
         result = pow(zgh / (zgh + delta), z - constants::half<T>());
      }
      // Split the calculation up to avoid spurious overflow:
      result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
   }
   result *= pow(constants::e<T>() / (zgh + delta), delta);
   return result;
}
//
// And again without Lanczos support this time:
//
template <class T, class Policy>
T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
{
   BOOST_MATH_STD_USING
   //
   // The upper gamma fraction is *very* slow for z < 6, actually it's very
   // slow to converge everywhere but recursing until z > 6 gets rid of the
   // worst of it's behaviour.
   //
   T prefix = 1;
   T zd = z + delta;
   while((zd < 6) && (z < 6))
   {
      prefix /= z;
      prefix *= zd;
      z += 1;
      zd += 1;
   }
   if(delta < 10)
   {
      prefix *= exp(-z * boost::math::log1p(delta / z, pol));
   }
   else
   {
      prefix *= pow(z / zd, z);
   }
   prefix *= pow(constants::e<T>() / zd, delta);
   T sum = detail::lower_gamma_series(z, z, pol) / z;
   sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
   T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
   sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
   sum /= sumd;
   if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
      return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
   return sum * prefix;
}

template <class T, class Policy>
T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
{
   BOOST_MATH_STD_USING

   if((z <= 0) || (z + delta <= 0))
   {
      // This isn't very sofisticated, or accurate, but it does work:
      return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
   }

   if(floor(delta) == delta)
   {
      if(floor(z) == z)
      {
         //
         // Both z and delta are integers, see if we can just use table lookup
         // of the factorials to get the result:
         //
         if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
         {
            return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
         }
      }
      if(fabs(delta) < 20)
      {
         //
         // delta is a small integer, we can use a finite product:
         //
         if(delta == 0)
            return 1;
         if(delta < 0)
         {
            z -= 1;
            T result = z;
            while(0 != (delta += 1))
            {
               z -= 1;
               result *= z;
            }
            return result;
         }
         else
         {
            T result = 1 / z;
            while(0 != (delta -= 1))
            {
               z += 1;
               result /= z;
            }
            return result;
         }
      }
   }
   typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
   return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
}

template <class T, class Policy>
T tgamma_ratio_imp(T x, T y, const Policy& pol)
{
   BOOST_MATH_STD_USING

   if((x <= 0) || (boost::math::isinf)(x))
      return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
   if((y <= 0) || (boost::math::isinf)(y))
      return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);

   if(x <= tools::min_value<T>())
   {
      // Special case for denorms...Ugh.
      T shift = ldexp(T(1), tools::digits<T>());
      return shift * tgamma_ratio_imp(T(x * shift), y, pol);
   }

   if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
   {
      // Rather than subtracting values, lets just call the gamma functions directly:
      return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
   }
   T prefix = 1;
   if(x < 1)
   {
      if(y < 2 * max_factorial<T>::value)
      {
         // We need to sidestep on x as well, otherwise we'll underflow
         // before we get to factor in the prefix term:
         prefix /= x;
         x += 1;
         while(y >=  max_factorial<T>::value)
         {
            y -= 1;
            prefix /= y;
         }
         return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
      }
      //
      // result is almost certainly going to underflow to zero, try logs just in case:
      //
      return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
   }
   if(y < 1)
   {
      if(x < 2 * max_factorial<T>::value)
      {
         // We need to sidestep on y as well, otherwise we'll overflow
         // before we get to factor in the prefix term:
         prefix *= y;
         y += 1;
         while(x >= max_factorial<T>::value)
         {
            x -= 1;
            prefix *= x;
         }
         return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
      }
      //
      // Result will almost certainly overflow, try logs just in case:
      //
      return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
   }
   //
   // Regular case, x and y both large and similar in magnitude:
   //
   return boost::math::tgamma_delta_ratio(x, y - x, pol);
}

template <class T, class Policy>
T gamma_p_derivative_imp(T a, T x, const Policy& pol)
{
   BOOST_MATH_STD_USING
   //
   // Usual error checks first:
   //
   if(a <= 0)
      return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
   if(x < 0)
      return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
   //
   // Now special cases:
   //
   if(x == 0)
   {
      return (a > 1) ? 0 :
         (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
   }
   //
   // Normal case:
   //
   typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
   T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
   if((x < 1) && (tools::max_value<T>() * x < f1))
   {
      // overflow:
      return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
   }
   if(f1 == 0)
   {
      // Underflow in calculation, use logs instead:
      f1 = a * log(x) - x - lgamma(a, pol) - log(x);
      f1 = exp(f1);
   }
   else
      f1 /= x;

   return f1;
}

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   tgamma(T z, const Policy& /* pol */, const mpl::true_)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;
   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)");
}

template <class T, class Policy>
struct igamma_initializer
{
   struct init
   {
      init()
      {
         typedef typename policies::precision<T, Policy>::type precision_type;

         typedef typename mpl::if_<
            mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
            mpl::greater<precision_type, mpl::int_<113> > >,
            mpl::int_<0>,
            typename mpl::if_<
               mpl::less_equal<precision_type, mpl::int_<53> >,
               mpl::int_<53>,
               typename mpl::if_<
                  mpl::less_equal<precision_type, mpl::int_<64> >,
                  mpl::int_<64>,
                  mpl::int_<113>
               >::type
            >::type
         >::type tag_type;

         do_init(tag_type());
      }
      template <int N>
      static void do_init(const mpl::int_<N>&)
      {
         // If std::numeric_limits<T>::digits is zero, we must not call
         // our inituialization code here as the precision presumably
         // varies at runtime, and will not have been set yet.  Plus the
         // code requiring initialization isn't called when digits == 0.
         if(std::numeric_limits<T>::digits)
         {
            boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
         }
      }
      static void do_init(const mpl::int_<53>&){}
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

template <class T, class Policy>
const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;

template <class T, class Policy>
struct lgamma_initializer
{
   struct init
   {
      init()
      {
         typedef typename policies::precision<T, Policy>::type precision_type;
         typedef typename mpl::if_<
            mpl::and_<
               mpl::less_equal<precision_type, mpl::int_<64> >,
               mpl::greater<precision_type, mpl::int_<0> >
            >,
            mpl::int_<64>,
            typename mpl::if_<
               mpl::and_<
                  mpl::less_equal<precision_type, mpl::int_<113> >,
                  mpl::greater<precision_type, mpl::int_<0> >
               >,
               mpl::int_<113>, mpl::int_<0> >::type
             >::type tag_type;
         do_init(tag_type());
      }
      static void do_init(const mpl::int_<64>&)
      {
         boost::math::lgamma(static_cast<T>(2.5), Policy());
         boost::math::lgamma(static_cast<T>(1.25), Policy());
         boost::math::lgamma(static_cast<T>(1.75), Policy());
      }
      static void do_init(const mpl::int_<113>&)
      {
         boost::math::lgamma(static_cast<T>(2.5), Policy());
         boost::math::lgamma(static_cast<T>(1.25), Policy());
         boost::math::lgamma(static_cast<T>(1.5), Policy());
         boost::math::lgamma(static_cast<T>(1.75), Policy());
      }
      static void do_init(const mpl::int_<0>&)
      {
      }
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

template <class T, class Policy>
const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;

template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   igamma_initializer<value_type, forwarding_policy>::force_instantiate();

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
      detail::gamma_incomplete_imp(static_cast<value_type>(a),
      static_cast<value_type>(z), false, true,
      forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
}

template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   tgamma(T1 a, T2 z, const mpl::false_ tag)
{
   return tgamma(a, z, policies::policy<>(), tag);
}

} // namespace detail

template <class T>
inline typename tools::promote_args<T>::type
   tgamma(T z)
{
   return tgamma(z, policies::policy<>());
}

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   lgamma(T z, int* sign, const Policy&)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");
}

template <class T>
inline typename tools::promote_args<T>::type
   lgamma(T z, int* sign)
{
   return lgamma(z, sign, policies::policy<>());
}

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   lgamma(T x, const Policy& pol)
{
   return ::boost::math::lgamma(x, 0, pol);
}

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

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   tgamma1pm1(T z, const Policy& /* pol */)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   return policies::checked_narrowing_cast<typename remove_cv<result_type>::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)");
}

template <class T>
inline typename tools::promote_args<T>::type
   tgamma1pm1(T z)
{
   return tgamma1pm1(z, policies::policy<>());
}

//
// Full upper incomplete gamma:
//
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   tgamma(T1 a, T2 z)
{
   //
   // Type T2 could be a policy object, or a value, select the
   // right overload based on T2:
   //
   typedef typename policies::is_policy<T2>::type maybe_policy;
   return detail::tgamma(a, z, maybe_policy());
}
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   tgamma(T1 a, T2 z, const Policy& pol)
{
   return detail::tgamma(a, z, pol, mpl::false_());
}
//
// Full lower incomplete gamma:
//
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   tgamma_lower(T1 a, T2 z, const Policy&)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
      detail::gamma_incomplete_imp(static_cast<value_type>(a),
      static_cast<value_type>(z), false, false,
      forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   tgamma_lower(T1 a, T2 z)
{
   return tgamma_lower(a, z, policies::policy<>());
}
//
// Regularised upper incomplete gamma:
//
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   gamma_q(T1 a, T2 z, const Policy& /* pol */)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
      detail::gamma_incomplete_imp(static_cast<value_type>(a),
      static_cast<value_type>(z), true, true,
      forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   gamma_q(T1 a, T2 z)
{
   return gamma_q(a, z, policies::policy<>());
}
//
// Regularised lower incomplete gamma:
//
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   gamma_p(T1 a, T2 z, const Policy&)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
      detail::gamma_incomplete_imp(static_cast<value_type>(a),
      static_cast<value_type>(z), true, false,
      forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   gamma_p(T1 a, T2 z)
{
   return gamma_p(a, z, policies::policy<>());
}

// ratios of gamma functions:
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, 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;

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   tgamma_delta_ratio(T1 z, T2 delta)
{
   return tgamma_delta_ratio(z, delta, policies::policy<>());
}
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   tgamma_ratio(T1 a, T2 b, const Policy&)
{
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, 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;

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   tgamma_ratio(T1 a, T2 b)
{
   return tgamma_ratio(a, b, policies::policy<>());
}

template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type
   gamma_p_derivative(T1 a, T2 x, const Policy&)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, 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;

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), forwarding_policy()), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
   gamma_p_derivative(T1 a, T2 x)
{
   return gamma_p_derivative(a, x, policies::policy<>());
}

} // namespace math
} // namespace boost

#ifdef BOOST_MSVC
# pragma warning(pop)
#endif

#include <boost/math/special_functions/detail/igamma_inverse.hpp>
#include <boost/math/special_functions/detail/gamma_inva.hpp>
#include <boost/math/special_functions/erf.hpp>

#endif // BOOST_MATH_SF_GAMMA_HPP


code:
:doctype: manpage
:man source: AsciiDoc
:man version: {VERSION}
:man manual: WSPRD Manual
= wsprd(1)

== NAME

wsprd - is a decoder for K1JT's Weak Signal Propagation Reporter (WSPR) mode.

== SYNOPSIS

*wsprd* ['OPTIONS']

== DESCRIPTION

*wsprd* - The program is written in C and is a command-line program that reads
from a .c2 file or .wav file and writes output to the console. It is used by
WSJT-X for wspr-mode decoding.

== OPTIONS
*-a <path>*:: Path to writeable data files, default="."

*-c <file>*:: Write .c2 file at the end of the first pass

*-e x*:: x is transceiver dial frequency error in Hz

*-f x*:: x is transceiver dial frequency in MHz

*-H* :: do not use, or update the hash table

*-m* :: decode wspr-15 .wav file

*-q* :: quick mode - does not dig deep for weak signals

*-s* :: single pass mode, no subtraction (same as original wsprd)

*-v* :: verbose mode, shows duplicate decodings

*-w* :: wideband mode - decode signals within {plus}/- 150 Hz of center

*-z x*:: x is fano metric table bias, default is 0.42

The Infile can be either .wav or .c2, for example:

-----
./wsprd -wf 14.0956 140709_2258.wav
-----

*NOTE* for .c2 files, the frequency within the file overrides the command
line value.

== FEATURES
* By default, *wsprd* reports signals that are within {plus}/- 110 Hz of the
subband center frequency. The wideband option (-w) extends this to {plus}/- 150 Hz.

* *wsprd* maintains a hashtable and will decode all three types of wspr
messages. An option (-H) is available to turn off use of the hashtable.

* The symbols are decoded using Phil Karn's sequential decoder routine, fano.c

== NOTES
. This program attempts to maximize the number of successful decodes per transmit
interval by trying to decode virtually every peak in the averaged spectrum.
The program also implements two-pass decoding, whereby signals that are successfully
decoded are subtracted one-by-one during the first decoding pass. Then, the
decoder is run again. In many cases the subtraction process will uncover signals
that can then be successfully decoded on the second pass.

. There will be occasional duplicate decodes when two closely spaced
peaks come from the same signal. The program removes dupes based on callsign
and frequency. Two decodes that have the same callsign and estimated frequencies
that are within 1 Hz will be treated as decodes of the same signal. This
dupechecking is turned off with the -v flag.

== AUTHORS

Joe Taylor, K1JT and Steven Franks, K9AN

== COPYING

*WSPRD* is Copyright (C) 2015 by Joseph H. Taylor, Jr., K1JT,
and Steven Franke, K9AN, with contributions from additional authors.
*WSPRD* is Open Source software, licensed under the GNU General Public
License (GPLv3).

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.


code:
#ifndef BOOST_ARCHIVE_DETAIL_AUTO_LINK_ARCHIVE_HPP
#define BOOST_ARCHIVE_DETAIL_AUTO_LINK_ARCHIVE_HPP

// MS compatible compilers support #pragma once
#if defined(_MSC_VER)
# pragma once
#endif

/////////1/////////2/////////3/////////4/////////5/////////6/////////7/////////8
//  auto_link_archive.hpp
//
//  (c) Copyright Robert Ramey 2004
//  Use, modification, and distribution is 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)

//  See library home page at http://www.boost.org/libs/serialization

//----------------------------------------------------------------------------//

// This header implements separate compilation features as described in
// http://www.boost.org/more/separate_compilation.html

//  enable automatic library variant selection  ------------------------------//

#include <boost/archive/detail/decl.hpp>

#if !defined(BOOST_ALL_NO_LIB) && !defined(BOOST_SERIALIZATION_NO_LIB) \
&&  !defined(BOOST_ARCHIVE_SOURCE) && !defined(BOOST_WARCHIVE_SOURCE)  \
&&  !defined(BOOST_SERIALIZATION_SOURCE)

    // Set the name of our library, this will get undef'ed by auto_link.hpp
    // once it's done with it:
    //
    #define BOOST_LIB_NAME boost_serialization
    //
    // If we're importing code from a dll, then tell auto_link.hpp about it:
    //
    #if defined(BOOST_ALL_DYN_LINK) || defined(BOOST_SERIALIZATION_DYN_LINK)
    #  define BOOST_DYN_LINK
    #endif
    //
    // And include the header that does the work:
    //
    #include <boost/config/auto_link.hpp>
#endif  // auto-linking disabled

#endif // BOOST_ARCHIVE_DETAIL_AUTO_LINK_ARCHIVE_HPP


code:
//////////////////////////////////////////////////////////////////////////////
//
// (C) Copyright David Abrahams, Vicente Botet 2009.
// (C) Copyright Ion Gaztanaga 2009-2012.
// 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
//! A general library header that includes
//! the rest of top-level headers.

#ifndef BOOST_MOVE_MOVE_HPP
#define BOOST_MOVE_MOVE_HPP

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

#include <boost/move/detail/config_begin.hpp>
#include <boost/move/utility.hpp>
#include <boost/move/iterator.hpp>
#include <boost/move/traits.hpp>
#include <boost/move/algorithm.hpp>
#include <boost/move/detail/config_end.hpp>

#endif //#ifndef BOOST_MOVE_MOVE_HPP


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_TOOLS_NEWTON_SOLVER_HPP
#define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <utility>
#include <boost/config/no_tr1/cmath.hpp>
#include <stdexcept>

#include <boost/math/tools/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/assert.hpp>
#include <boost/throw_exception.hpp>

#ifdef BOOST_MSVC
#pragma warning(push)
#pragma warning(disable: 4512)
#endif
#include <boost/math/tools/tuple.hpp>
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif

#include <boost/math/special_functions/sign.hpp>
#include <boost/math/tools/toms748_solve.hpp>
#include <boost/math/policies/error_handling.hpp>

namespace boost{ namespace math{ namespace tools{

namespace detail{

namespace dummy{

   template<int n, class T>
   typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T);
}

template <class Tuple, class T>
void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T)
{
   using dummy::get;
   // Use ADL to find the right overload for get:
   a = get<0>(t);
   b = get<1>(t);
}
template <class Tuple, class T>
void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T)
{
   using dummy::get;
   // Use ADL to find the right overload for get:
   a = get<0>(t);
   b = get<1>(t);
   c = get<2>(t);
}

template <class Tuple, class T>
inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
{
   using dummy::get;
   // Rely on ADL to find the correct overload of get:
   val = get<0>(t);
}

template <class T, class U, class V>
inline void unpack_tuple(const std::pair<T, U>& p, V& a, V& b) BOOST_MATH_NOEXCEPT(T)
{
   a = p.first;
   b = p.second;
}
template <class T, class U, class V>
inline void unpack_0(const std::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T)
{
   a = p.first;
}

template <class F, class T>
void handle_zero_derivative(F f,
                            T& last_f0,
                            const T& f0,
                            T& delta,
                            T& result,
                            T& guess,
                            const T& min,
                            const T& max) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   if(last_f0 == 0)
   {
      // this must be the first iteration, pretend that we had a
      // previous one at either min or max:
      if(result == min)
      {
         guess = max;
      }
      else
      {
         guess = min;
      }
      unpack_0(f(guess), last_f0);
      delta = guess - result;
   }
   if(sign(last_f0) * sign(f0) < 0)
   {
      // we've crossed over so move in opposite direction to last step:
      if(delta < 0)
      {
         delta = (result - min) / 2;
      }
      else
      {
         delta = (result - max) / 2;
      }
   }
   else
   {
      // move in same direction as last step:
      if(delta < 0)
      {
         delta = (result - max) / 2;
      }
      else
      {
         delta = (result - min) / 2;
      }
   }
}

} // namespace

template <class F, class T, class Tol, class Policy>
std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<Policy>::value && BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   T fmin = f(min);
   T fmax = f(max);
   if(fmin == 0)
   {
      max_iter = 2;
      return std::make_pair(min, min);
   }
   if(fmax == 0)
   {
      max_iter = 2;
      return std::make_pair(max, max);
   }

   //
   // Error checking:
   //
   static const char* function = "boost::math::tools::bisect<%1%>";
   if(min >= max)
   {
      return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
         "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol));
   }
   if(fmin * fmax >= 0)
   {
      return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
         "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = %1%).", fmin, pol));
   }

   //
   // Three function invocations so far:
   //
   boost::uintmax_t count = max_iter;
   if(count < 3)
      count = 0;
   else
      count -= 3;

   while(count && (0 == tol(min, max)))
   {
      T mid = (min + max) / 2;
      T fmid = f(mid);
      if((mid == max) || (mid == min))
         break;
      if(fmid == 0)
      {
         min = max = mid;
         break;
      }
      else if(sign(fmid) * sign(fmin) < 0)
      {
         max = mid;
         fmax = fmid;
      }
      else
      {
         min = mid;
         fmin = fmid;
      }
      --count;
   }

   max_iter -= count;

#ifdef BOOST_MATH_INSTRUMENT
   std::cout << "Bisection iteration, final count = " << max_iter << std::endl;

   static boost::uintmax_t max_count = 0;
   if(max_iter > max_count)
   {
      max_count = max_iter;
      std::cout << "Maximum iterations: " << max_iter << std::endl;
   }
#endif

   return std::make_pair(min, max);
}

template <class F, class T, class Tol>
inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter)  BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   return bisect(f, min, max, tol, max_iter, policies::policy<>());
}

template <class F, class T, class Tol>
inline std::pair<T, T> bisect(F f, T min, T max, Tol tol) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
   return bisect(f, min, max, tol, m, policies::policy<>());
}

template <class F, class T>
T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   BOOST_MATH_STD_USING

   T f0(0), f1, last_f0(0);
   T result = guess;

   T factor = static_cast<T>(ldexp(1.0, 1 - digits));
   T delta = tools::max_value<T>();
   T delta1 = tools::max_value<T>();
   T delta2 = tools::max_value<T>();

   boost::uintmax_t count(max_iter);

   do{
      last_f0 = f0;
      delta2 = delta1;
      delta1 = delta;
      detail::unpack_tuple(f(result), f0, f1);
      --count;
      if(0 == f0)
         break;
      if(f1 == 0)
      {
         // Oops zero derivative!!!
#ifdef BOOST_MATH_INSTRUMENT
         std::cout << "Newton iteration, zero derivative found" << std::endl;
#endif
         detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
      }
      else
      {
         delta = f0 / f1;
      }
#ifdef BOOST_MATH_INSTRUMENT
      std::cout << "Newton iteration, delta = " << delta << std::endl;
#endif
      if(fabs(delta * 2) > fabs(delta2))
      {
         // last two steps haven't converged, try bisection:
         delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
      }
      guess = result;
      result -= delta;
      if(result <= min)
      {
         delta = 0.5F * (guess - min);
         result = guess - delta;
         if((result == min) || (result == max))
            break;
      }
      else if(result >= max)
      {
         delta = 0.5F * (guess - max);
         result = guess - delta;
         if((result == min) || (result == max))
            break;
      }
      // update brackets:
      if(delta > 0)
         max = guess;
      else
         min = guess;
   }while(count && (fabs(result * factor) < fabs(delta)));

   max_iter -= count;

#ifdef BOOST_MATH_INSTRUMENT
   std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;

   static boost::uintmax_t max_count = 0;
   if(max_iter > max_count)
   {
      max_count = max_iter;
      std::cout << "Maximum iterations: " << max_iter << std::endl;
   }
#endif

   return result;
}

template <class F, class T>
inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
   return newton_raphson_iterate(f, guess, min, max, digits, m);
}

namespace detail{

   struct halley_step
   {
      template <class T>
      static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
      {
         using std::fabs;
         T denom = 2 * f0;
         T num = 2 * f1 - f0 * (f2 / f1);
         T delta;

         BOOST_MATH_INSTRUMENT_VARIABLE(denom);
         BOOST_MATH_INSTRUMENT_VARIABLE(num);

         if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
         {
            // possible overflow, use Newton step:
            delta = f0 / f1;
         }
         else
            delta = denom / num;
         return delta;
      }
   };

   template <class Stepper, class F, class T>
   T second_order_root_finder(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
   {
      BOOST_MATH_STD_USING

         T f0(0), f1, f2;
      T result = guess;

      T factor = static_cast<T>(ldexp(1.0, 1 - digits));
      T delta = (std::max)(T(10000000 * guess), T(10000000));  // arbitarily large delta
      T last_f0 = 0;
      T delta1 = delta;
      T delta2 = delta;

      bool out_of_bounds_sentry = false;

#ifdef BOOST_MATH_INSTRUMENT
      std::cout << "Second order root iteration, limit = " << factor << std::endl;
#endif

      boost::uintmax_t count(max_iter);

      do{
         last_f0 = f0;
         delta2 = delta1;
         delta1 = delta;
         detail::unpack_tuple(f(result), f0, f1, f2);
         --count;

         BOOST_MATH_INSTRUMENT_VARIABLE(f0);
         BOOST_MATH_INSTRUMENT_VARIABLE(f1);
         BOOST_MATH_INSTRUMENT_VARIABLE(f2);

         if(0 == f0)
            break;
         if(f1 == 0)
         {
            // Oops zero derivative!!!
#ifdef BOOST_MATH_INSTRUMENT
            std::cout << "Second order root iteration, zero derivative found" << std::endl;
#endif
            detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
         }
         else
         {
            if(f2 != 0)
            {
               delta = Stepper::step(result, f0, f1, f2);
               if(delta * f1 / f0 < 0)
               {
                  // Oh dear, we have a problem as Newton and Halley steps
                  // disagree about which way we should move.  Probably
                  // there is cancelation error in the calculation of the
                  // Halley step, or else the derivatives are so small
                  // that their values are basically trash.  We will move
                  // in the direction indicated by a Newton step, but
                  // by no more than twice the current guess value, otherwise
                  // we can jump way out of bounds if we're not careful.
                  // See https://svn.boost.org/trac/boost/ticket/8314.
                  delta = f0 / f1;
                  if(fabs(delta) > 2 * fabs(guess))
                     delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
               }
            }
            else
               delta = f0 / f1;
         }
#ifdef BOOST_MATH_INSTRUMENT
         std::cout << "Second order root iteration, delta = " << delta << std::endl;
#endif
         T convergence = fabs(delta / delta2);
         if((convergence > 0.8) && (convergence < 2))
         {
            // last two steps haven't converged, try bisection:
            delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
            if(fabs(delta) > result)
               delta = sign(delta) * result; // protect against huge jumps!
            // reset delta2 so that this branch will *not* be taken on the
            // next iteration:
            delta2 = delta * 3;
            BOOST_MATH_INSTRUMENT_VARIABLE(delta);
         }
         guess = result;
         result -= delta;
         BOOST_MATH_INSTRUMENT_VARIABLE(result);

         // check for out of bounds step:
         if(result < min)
         {
            T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
            if(fabs(diff) < 1)
               diff = 1 / diff;
            if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
            {
               // Only a small out of bounds step, lets assume that the result
               // is probably approximately at min:
               delta = 0.99f * (guess - min);
               result = guess - delta;
               out_of_bounds_sentry = true; // only take this branch once!
            }
            else
            {
               delta = (guess - min) / 2;
               result = guess - delta;
               if((result == min) || (result == max))
                  break;
            }
         }
         else if(result > max)
         {
            T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
            if(fabs(diff) < 1)
               diff = 1 / diff;
            if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
            {
               // Only a small out of bounds step, lets assume that the result
               // is probably approximately at min:
               delta = 0.99f * (guess - max);
               result = guess - delta;
               out_of_bounds_sentry = true; // only take this branch once!
            }
            else
            {
               delta = (guess - max) / 2;
               result = guess - delta;
               if((result == min) || (result == max))
                  break;
            }
         }
         // update brackets:
         if(delta > 0)
            max = guess;
         else
            min = guess;
      } while(count && (fabs(result * factor) < fabs(delta)));

      max_iter -= count;

#ifdef BOOST_MATH_INSTRUMENT
      std::cout << "Second order root iteration, final count = " << max_iter << std::endl;
#endif

      return result;
   }

}

template <class F, class T>
T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter);
}

template <class F, class T>
inline T halley_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
   return halley_iterate(f, guess, min, max, digits, m);
}

namespace detail{

   struct schroder_stepper
   {
      template <class T>
      static T step(const T& x, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
      {
         T ratio = f0 / f1;
         T delta;
         if(ratio / x < 0.1)
         {
            delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
            // check second derivative doesn't over compensate:
            if(delta * ratio < 0)
               delta = ratio;
         }
         else
            delta = ratio;  // fall back to Newton iteration.
         return delta;
      }
   };

}

template <class F, class T>
T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
}

template <class F, class T>
inline T schroder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
   return schroder_iterate(f, guess, min, max, digits, m);
}
//
// These two are the old spelling of this function, retained for backwards compatibity just in case:
//
template <class F, class T>
T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
}

template <class F, class T>
inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
   boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
   return schroder_iterate(f, guess, min, max, digits, m);
}

} // namespace tools
} // namespace math
} // namespace boost

#endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP


code:
#ifndef BOOST_MPL_LIST_LIST30_HPP_INCLUDED
#define BOOST_MPL_LIST_LIST30_HPP_INCLUDED

// Copyright Aleksey Gurtovoy 2000-2004
//
// 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$

#if !defined(BOOST_MPL_PREPROCESSING_MODE)
#   include <boost/mpl/list/list20.hpp>
#endif

#include <boost/mpl/aux_/config/use_preprocessed.hpp>

#if !defined(BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS) \
 && !defined(BOOST_MPL_PREPROCESSING_MODE)

#   define BOOST_MPL_PREPROCESSED_HEADER list30.hpp
#   include <boost/mpl/list/aux_/include_preprocessed.hpp>

#else

#   include <boost/preprocessor/iterate.hpp>

namespace boost { namespace mpl {

#   define BOOST_PP_ITERATION_PARAMS_1 \
    (3,(21, 30, <boost/mpl/list/aux_/numbered.hpp>))
#   include BOOST_PP_ITERATE()

}}

#endif // BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS

#endif // BOOST_MPL_LIST_LIST30_HPP_INCLUDED


code:
#ifndef BOOST_MPL_AUX_NUMERIC_CAST_HPP_INCLUDED
#define BOOST_MPL_AUX_NUMERIC_CAST_HPP_INCLUDED

// Copyright Aleksey Gurtovoy 2003-2004
//
// 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$

#include <boost/mpl/numeric_cast.hpp>
#include <boost/mpl/apply_wrap.hpp>
#include <boost/mpl/aux_/config/forwarding.hpp>

namespace boost { namespace mpl { namespace aux {

template<
      typename F
    , typename Tag1
    , typename Tag2
    >
struct cast1st_impl
{
    template< typename N1, typename N2 > struct apply
#if !defined(BOOST_MPL_CFG_NO_NESTED_FORWARDING)
        : apply_wrap2<
              F
            , typename apply_wrap1< BOOST_MPL_AUX_NUMERIC_CAST<Tag1,Tag2>,N1 >::type
            , N2
            >
    {
#else
    {
    typedef typename apply_wrap2<
              F
            , typename apply_wrap1< BOOST_MPL_AUX_NUMERIC_CAST<Tag1,Tag2>,N1 >::type
            , N2
            >::type type;
#endif
    };
};

template<
      typename F
    , typename Tag1
    , typename Tag2
    >
struct cast2nd_impl
{
    template< typename N1, typename N2 > struct apply
#if !defined(BOOST_MPL_CFG_NO_NESTED_FORWARDING)
        : apply_wrap2<
              F
            , N1
            , typename apply_wrap1< BOOST_MPL_AUX_NUMERIC_CAST<Tag2,Tag1>,N2 >::type
            >
    {
#else
    {
        typedef typename apply_wrap2<
              F
            , N1
            , typename apply_wrap1< BOOST_MPL_AUX_NUMERIC_CAST<Tag2,Tag1>,N2 >::type
            >::type type;
#endif
    };
};

}}}

#endif // BOOST_MPL_AUX_NUMERIC_CAST_HPP_INCLUDED


code:
// Copyright Vladimir Prus 2004.
// 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_PROGRAM_OPTIONS_POSITIONAL_OPTIONS_VP_2004_03_02
#define BOOST_PROGRAM_OPTIONS_POSITIONAL_OPTIONS_VP_2004_03_02

#include <boost/program_options/config.hpp>

#include <vector>
#include <string>

#if defined(BOOST_MSVC)
#   pragma warning (push)
#   pragma warning (disable:4251) // class 'std::vector<_Ty>' needs to have dll-interface to be used by clients of class 'boost::program_options::positional_options_description'
#endif

namespace boost { namespace program_options {

    /** Describes positional options.

        The class allows to guess option names for positional options, which
        are specified on the command line and are identified by the position.
        The class uses the information provided by the user to associate a name
        with every positional option, or tell that no name is known.

        The primary assumption is that only the relative order of the
        positional options themselves matters, and that any interleaving
        ordinary options don't affect interpretation of positional options.

        The user initializes the class by specifying that first N positional
        options should be given the name X1, following M options should be given
        the name X2 and so on.
    */
    class BOOST_PROGRAM_OPTIONS_DECL positional_options_description {
    public:
        positional_options_description();

        /** Species that up to 'max_count' next positional options
            should be given the 'name'. The value of '-1' means 'unlimited'.
            No calls to 'add' can be made after call with 'max_value' equal to
            '-1'.
        */
        positional_options_description&
        add(const char* name, int max_count);

        /** Returns the maximum number of positional options that can
            be present. Can return (numeric_limits<unsigned>::max)() to
            indicate unlimited number. */
        unsigned max_total_count() const;

        /** Returns the name that should be associated with positional
            options at 'position'.
            Precondition: position < max_total_count()
        */
        const std::string& name_for_position(unsigned position) const;

    private:
        // List of names corresponding to the positions. If the number of
        // positions is unlimited, then the last name is stored in
        // m_trailing;
        std::vector<std::string> m_names;
        std::string m_trailing;
    };

}}

#if defined(BOOST_MSVC)
#   pragma warning (pop)
#endif

#endif