code:
//  Copyright John Maddock 2007.
//  Copyright Paul A. Bristow 2007
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_SF_DETAIL_INV_T_HPP
#define BOOST_MATH_SF_DETAIL_INV_T_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/cbrt.hpp>
#include <boost/math/special_functions/round.hpp>
#include <boost/math/special_functions/trunc.hpp>

namespace boost{ namespace math{ namespace detail{

//
// The main method used is due to Hill:
//
// G. W. Hill, Algorithm 396, Student's t-Quantiles,
// Communications of the ACM, 13(10): 619-620, Oct., 1970.
//
template <class T, class Policy>
T inverse_students_t_hill(T ndf, T u, const Policy& pol)
{
   BOOST_MATH_STD_USING
   BOOST_ASSERT(u <= 0.5);

   T a, b, c, d, q, x, y;

   if (ndf > 1e20f)
      return -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();

   a = 1 / (ndf - 0.5f);
   b = 48 / (a * a);
   c = ((20700 * a / b - 98) * a - 16) * a + 96.36f;
   d = ((94.5f / (b + c) - 3) / b + 1) * sqrt(a * constants::pi<T>() / 2) * ndf;
   y = pow(d * 2 * u, 2 / ndf);

   if (y > (0.05f + a))
   {
      //
      // Asymptotic inverse expansion about normal:
      //
      x = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
      y = x * x;

      if (ndf < 5)
         c += 0.3f * (ndf - 4.5f) * (x + 0.6f);
      c += (((0.05f * d * x - 5) * x - 7) * x - 2) * x + b;
      y = (((((0.4f * y + 6.3f) * y + 36) * y + 94.5f) / c - y - 3) / b + 1) * x;
      y = boost::math::expm1(a * y * y, pol);
   }
   else
   {
      y = static_cast<T>(((1 / (((ndf + 6) / (ndf * y) - 0.089f * d - 0.822f)
              * (ndf + 2) * 3) + 0.5 / (ndf + 4)) * y - 1)
              * (ndf + 1) / (ndf + 2) + 1 / y);
   }
   q = sqrt(ndf * y);

   return -q;
}
//
// Tail and body series are due to Shaw:
//
// www.mth.kcl.ac.uk/~shaww/web_page/papers/Tdistribution06.pdf
//
// Shaw, W.T., 2006, "Sampling Student's T distribution - use of
// the inverse cumulative distribution function."
// Journal of Computational Finance, Vol 9 Issue 4, pp 37-73, Summer 2006
//
template <class T, class Policy>
T inverse_students_t_tail_series(T df, T v, const Policy& pol)
{
   BOOST_MATH_STD_USING
   // Tail series expansion, see section 6 of Shaw's paper.
   // w is calculated using Eq 60:
   T w = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
      * sqrt(df * constants::pi<T>()) * v;
   // define some variables:
   T np2 = df + 2;
   T np4 = df + 4;
   T np6 = df + 6;
   //
   // Calculate the coefficients d(k), these depend only on the
   // number of degrees of freedom df, so at least in theory
   // we could tabulate these for fixed df, see p15 of Shaw:
   //
   T d[7] = { 1, };
   d[1] = -(df + 1) / (2 * np2);
   np2 *= (df + 2);
   d[2] = -df * (df + 1) * (df + 3) / (8 * np2 * np4);
   np2 *= df + 2;
   d[3] = -df * (df + 1) * (df + 5) * (((3 * df) + 7) * df -2) / (48 * np2 * np4 * np6);
   np2 *= (df + 2);
   np4 *= (df + 4);
   d[4] = -df * (df + 1) * (df + 7) *
      ( (((((15 * df) + 154) * df + 465) * df + 286) * df - 336) * df + 64 )
      / (384 * np2 * np4 * np6 * (df + 8));
   np2 *= (df + 2);
   d[5] = -df * (df + 1) * (df + 3) * (df + 9)
            * (((((((35 * df + 452) * df + 1573) * df + 600) * df - 2020) * df) + 928) * df -128)
            / (1280 * np2 * np4 * np6 * (df + 8) * (df + 10));
   np2 *= (df + 2);
   np4 *= (df + 4);
   np6 *= (df + 6);
   d[6] = -df * (df + 1) * (df + 11)
            * ((((((((((((945 * df) + 31506) * df + 425858) * df + 2980236) * df + 11266745) * df + 20675018) * df + 7747124) * df - 22574632) * df - 8565600) * df + 18108416) * df - 7099392) * df + 884736)
            / (46080 * np2 * np4 * np6 * (df + 8) * (df + 10) * (df +12));
   //
   // Now bring everthing together to provide the result,
   // this is Eq 62 of Shaw:
   //
   T rn = sqrt(df);
   T div = pow(rn * w, 1 / df);
   T power = div * div;
   T result = tools::evaluate_polynomial<7, T, T>(d, power);
   result *= rn;
   result /= div;
   return -result;
}

template <class T, class Policy>
T inverse_students_t_body_series(T df, T u, const Policy& pol)
{
   BOOST_MATH_STD_USING
   //
   // Body series for small N:
   //
   // Start with Eq 56 of Shaw:
   //
   T v = boost::math::tgamma_delta_ratio(df / 2, constants::half<T>(), pol)
      * sqrt(df * constants::pi<T>()) * (u - constants::half<T>());
   //
   // Workspace for the polynomial coefficients:
   //
   T c[11] = { 0, 1, };
   //
   // Figure out what the coefficients are, note these depend
   // only on the degrees of freedom (Eq 57 of Shaw):
   //
   T in = 1 / df;
   c[2] = static_cast<T>(0.16666666666666666667 + 0.16666666666666666667 * in);
   c[3] = static_cast<T>((0.0083333333333333333333 * in
      + 0.066666666666666666667) * in
      + 0.058333333333333333333);
   c[4] = static_cast<T>(((0.00019841269841269841270 * in
      + 0.0017857142857142857143) * in
      + 0.026785714285714285714) * in
      + 0.025198412698412698413);
   c[5] = static_cast<T>((((2.7557319223985890653e-6 * in
      + 0.00037477954144620811287) * in
      - 0.0011078042328042328042) * in
      + 0.010559964726631393298) * in
      + 0.012039792768959435626);
   c[6] = static_cast<T>(((((2.5052108385441718775e-8 * in
      - 0.000062705427288760622094) * in
      + 0.00059458674042007375341) * in
      - 0.0016095979637646304313) * in
      + 0.0061039211560044893378) * in
      + 0.0038370059724226390893);
   c[7] = static_cast<T>((((((1.6059043836821614599e-10 * in
      + 0.000015401265401265401265) * in
      - 0.00016376804137220803887) * in
      + 0.00069084207973096861986) * in
      - 0.0012579159844784844785) * in
      + 0.0010898206731540064873) * in
      + 0.0032177478835464946576);
   c[8] = static_cast<T>(((((((7.6471637318198164759e-13 * in
      - 3.9851014346715404916e-6) * in
      + 0.000049255746366361445727) * in
      - 0.00024947258047043099953) * in
      + 0.00064513046951456342991) * in
      - 0.00076245135440323932387) * in
      + 0.000033530976880017885309) * in
      + 0.0017438262298340009980);
   c[9] = static_cast<T>((((((((2.8114572543455207632e-15 * in
      + 1.0914179173496789432e-6) * in
      - 0.000015303004486655377567) * in
      + 0.000090867107935219902229) * in
      - 0.00029133414466938067350) * in
      + 0.00051406605788341121363) * in
      - 0.00036307660358786885787) * in
      - 0.00031101086326318780412) * in
      + 0.00096472747321388644237);
   c[10] = static_cast<T>(((((((((8.2206352466243297170e-18 * in
      - 3.1239569599829868045e-7) * in
      + 4.8903045291975346210e-6) * in
      - 0.000033202652391372058698) * in
      + 0.00012645437628698076975) * in
      - 0.00028690924218514613987) * in
      + 0.00035764655430568632777) * in
      - 0.00010230378073700412687) * in
      - 0.00036942667800009661203) * in
      + 0.00054229262813129686486);
   //
   // The result is then a polynomial in v (see Eq 56 of Shaw):
   //
   return tools::evaluate_odd_polynomial<11, T, T>(c, v);
}

template <class T, class Policy>
T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
{
   //
   // df = number of degrees of freedom.
   // u = probablity.
   // v = 1 - u.
   // l = lanczos type to use.
   //
   BOOST_MATH_STD_USING
   bool invert = false;
   T result = 0;
   if(pexact)
      *pexact = false;
   if(u > v)
   {
      // function is symmetric, invert it:
      std::swap(u, v);
      invert = true;
   }
   if((floor(df) == df) && (df < 20))
   {
      //
      // we have integer degrees of freedom, try for the special
      // cases first:
      //
      T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);

      switch(itrunc(df, Policy()))
      {
      case 1:
         {
            //
            // df = 1 is the same as the Cauchy distribution, see
            // Shaw Eq 35:
            //
            if(u == 0.5)
               result = 0;
            else
               result = -cos(constants::pi<T>() * u) / sin(constants::pi<T>() * u);
            if(pexact)
               *pexact = true;
            break;
         }
      case 2:
         {
            //
            // df = 2 has an exact result, see Shaw Eq 36:
            //
            result =(2 * u - 1) / sqrt(2 * u * v);
            if(pexact)
               *pexact = true;
            break;
         }
      case 4:
         {
            //
            // df = 4 has an exact result, see Shaw Eq 38 & 39:
            //
            T alpha = 4 * u * v;
            T root_alpha = sqrt(alpha);
            T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
            T x = sqrt(r - 4);
            result = u - 0.5f < 0 ? (T)-x : x;
            if(pexact)
               *pexact = true;
            break;
         }
      case 6:
         {
            //
            // We get numeric overflow in this area:
            //
            if(u < 1e-150)
               return (invert ? -1 : 1) * inverse_students_t_hill(df, u, pol);
            //
            // Newton-Raphson iteration of a polynomial case,
            // choice of seed value is taken from Shaw's online
            // supplement:
            //
            T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
            T b = boost::math::cbrt(a);
            static const T c = static_cast<T>(0.85498797333834849467655443627193);
            T p = 6 * (1 + c * (1 / b - 1));
            T p0;
            do{
               T p2 = p * p;
               T p4 = p2 * p2;
               T p5 = p * p4;
               p0 = p;
               // next term is given by Eq 41:
               p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
            }while(fabs((p - p0) / p) > tolerance);
            //
            // Use Eq 45 to extract the result:
            //
            p = sqrt(p - df);
            result = (u - 0.5f) < 0 ? (T)-p : p;
            break;
         }
#if 0
         //
         // These are Shaw's "exact" but iterative solutions
         // for even df, the numerical accuracy of these is
         // rather less than Hill's method, so these are disabled
         // for now, which is a shame because they are reasonably
         // quick to evaluate...
         //
      case 8:
         {
            //
            // Newton-Raphson iteration of a polynomial case,
            // choice of seed value is taken from Shaw's online
            // supplement:
            //
            static const T c8 = 0.85994765706259820318168359251872L;
            T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
            T b = pow(a, T(1) / 4);
            T p = 8 * (1 + c8 * (1 / b - 1));
            T p0 = p;
            do{
               T p5 = p * p;
               p5 *= p5 * p;
               p0 = p;
               // Next term is given by Eq 42:
               p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
            }while(fabs((p - p0) / p) > tolerance);
            //
            // Use Eq 45 to extract the result:
            //
            p = sqrt(p - df);
            result = (u - 0.5f) < 0 ? -p : p;
            break;
         }
      case 10:
         {
            //
            // Newton-Raphson iteration of a polynomial case,
            // choice of seed value is taken from Shaw's online
            // supplement:
            //
            static const T c10 = 0.86781292867813396759105692122285L;
            T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
            T b = pow(a, T(1) / 5);
            T p = 10 * (1 + c10 * (1 / b - 1));
            T p0;
            do{
               T p6 = p * p;
               p6 *= p6 * p6;
               p0 = p;
               // Next term given by Eq 43:
               p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) /
                  (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
            }while(fabs((p - p0) / p) > tolerance);
            //
            // Use Eq 45 to extract the result:
            //
            p = sqrt(p - df);
            result = (u - 0.5f) < 0 ? -p : p;
            break;
         }
#endif
      default:
         goto calculate_real;
      }
   }
   else
   {
calculate_real:
      if(df > 0x10000000)
      {
         result = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
         if((pexact) && (df >= 1e20))
            *pexact = true;
      }
      else if(df < 3)
      {
         //
         // Use a roughly linear scheme to choose between Shaw's
         // tail series and body series:
         //
         T crossover = 0.2742f - df * 0.0242143f;
         if(u > crossover)
         {
            result = boost::math::detail::inverse_students_t_body_series(df, u, pol);
         }
         else
         {
            result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
         }
      }
      else
      {
         //
         // Use Hill's method except in the exteme tails
         // where we use Shaw's tail series.
         // The crossover point is roughly exponential in -df:
         //
         T crossover = ldexp(1.0f, iround(T(df / -0.654f), typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
         if(u > crossover)
         {
            result = boost::math::detail::inverse_students_t_hill(df, u, pol);
         }
         else
         {
            result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
         }
      }
   }
   return invert ? (T)-result : result;
}

template <class T, class Policy>
inline T find_ibeta_inv_from_t_dist(T a, T p, T /*q*/, T* py, const Policy& pol)
{
   T u = p / 2;
   T v = 1 - u;
   T df = a * 2;
   T t = boost::math::detail::inverse_students_t(df, u, v, pol);
   *py = t * t / (df + t * t);
   return df / (df + t * t);
}

template <class T, class Policy>
inline T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::false_*)
{
   BOOST_MATH_STD_USING
   //
   // Need to use inverse incomplete beta to get
   // required precision so not so fast:
   //
   T probability = (p > 0.5) ? 1 - p : p;
   T t, x, y(0);
   x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y, pol);
   if(df * y > tools::max_value<T>() * x)
      t = policies::raise_overflow_error<T>("boost::math::students_t_quantile<%1%>(%1%,%1%)", 0, pol);
   else
      t = sqrt(df * y / x);
   //
   // Figure out sign based on the size of p:
   //
   if(p < 0.5)
      t = -t;
   return t;
}

template <class T, class Policy>
T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::true_*)
{
   BOOST_MATH_STD_USING
   bool invert = false;
   if((df < 2) && (floor(df) != df))
      return boost::math::detail::fast_students_t_quantile_imp(df, p, pol, static_cast<mpl::false_*>(0));
   if(p > 0.5)
   {
      p = 1 - p;
      invert = true;
   }
   //
   // Get an estimate of the result:
   //
   bool exact;
   T t = inverse_students_t(df, p, T(1-p), pol, &exact);
   if((t == 0) || exact)
      return invert ? -t : t; // can't do better!
   //
   // Change variables to inverse incomplete beta:
   //
   T t2 = t * t;
   T xb = df / (df + t2);
   T y = t2 / (df + t2);
   T a = df / 2;
   //
   // t can be so large that x underflows,
   // just return our estimate in that case:
   //
   if(xb == 0)
      return t;
   //
   // Get incomplete beta and it's derivative:
   //
   T f1;
   T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, pol, false, true, &f1)
      : ibeta_imp(constants::half<T>(), a, y, pol, true, true, &f1);

   // Get cdf from incomplete beta result:
   T p0 = f0 / 2  - p;
   // Get pdf from derivative:
   T p1 = f1 * sqrt(y * xb * xb * xb / df);
   //
   // Second derivative divided by p1:
   //
   // yacas gives:
   //
   // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
   //
   //  |                        | v + 1     |     |
   //  |                       -| ----- + 1 |     |
   //  |                        |   2       |     |
   // -|             |  2     |                   |
   //  |             | t      |                   |
   //  |             | -- + 1 |                   |
   //  | ( v + 1 ) * | v      |               * t |
   // ---------------------------------------------
   //                       v
   //
   // Which after some manipulation is:
   //
   // -p1 * t * (df + 1) / (t^2 + df)
   //
   T p2 = t * (df + 1) / (t * t + df);
   // Halley step:
   t = fabs(t);
   t += p0 / (p1 + p0 * p2 / 2);
   return !invert ? -t : t;
}

template <class T, class Policy>
inline T fast_students_t_quantile(T df, T p, const Policy& pol)
{
   typedef typename policies::evaluation<T, Policy>::type value_type;
   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   typedef mpl::bool_<
      (std::numeric_limits<T>::digits <= 53)
       &&
      (std::numeric_limits<T>::is_specialized)
       &&
      (std::numeric_limits<T>::radix == 2)
   > tag_type;
   return policies::checked_narrowing_cast<T, forwarding_policy>(fast_students_t_quantile_imp(static_cast<value_type>(df), static_cast<value_type>(p), pol, static_cast<tag_type*>(0)), "boost::math::students_t_quantile<%1%>(%1%,%1%,%1%)");
}

}}} // namespaces

#endif // BOOST_MATH_SF_DETAIL_INV_T_HPP


code:
//  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_SP_UC_FACTORIALS_HPP
#define BOOST_MATH_SP_UC_FACTORIALS_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/array.hpp>
#ifdef BOOST_MSVC
#pragma warning(push) // Temporary until lexical cast fixed.
#pragma warning(disable: 4127 4701)
#endif
#ifndef BOOST_MATH_NO_LEXICAL_CAST
#include <boost/lexical_cast.hpp>
#endif
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
#include <boost/config/no_tr1/cmath.hpp>
#include <boost/math/special_functions/math_fwd.hpp>

namespace boost { namespace math
{
// Forward declarations:
template <class T>
struct max_factorial;

// Definitions:
template <>
inline float unchecked_factorial<float>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(float))
{
   static const boost::array<float, 35> factorials = {{
      1.0F,
      1.0F,
      2.0F,
      6.0F,
      24.0F,
      120.0F,
      720.0F,
      5040.0F,
      40320.0F,
      362880.0F,
      3628800.0F,
      39916800.0F,
      479001600.0F,
      6227020800.0F,
      87178291200.0F,
      1307674368000.0F,
      20922789888000.0F,
      355687428096000.0F,
      6402373705728000.0F,
      121645100408832000.0F,
      0.243290200817664e19F,
      0.5109094217170944e20F,
      0.112400072777760768e22F,
      0.2585201673888497664e23F,
      0.62044840173323943936e24F,
      0.15511210043330985984e26F,
      0.403291461126605635584e27F,
      0.10888869450418352160768e29F,
      0.304888344611713860501504e30F,
      0.8841761993739701954543616e31F,
      0.26525285981219105863630848e33F,
      0.822283865417792281772556288e34F,
      0.26313083693369353016721801216e36F,
      0.868331761881188649551819440128e37F,
      0.29523279903960414084761860964352e39F,
   }};

   return factorials[i];
}

template <>
struct max_factorial<float>
{
   BOOST_STATIC_CONSTANT(unsigned, value = 34);
};

template <>
inline long double unchecked_factorial<long double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(long double))
{
   static const boost::array<long double, 171> factorials = {{
      1L,
      1L,
      2L,
      6L,
      24L,
      120L,
      720L,
      5040L,
      40320L,
      362880.0L,
      3628800.0L,
      39916800.0L,
      479001600.0L,
      6227020800.0L,
      87178291200.0L,
      1307674368000.0L,
      20922789888000.0L,
      355687428096000.0L,
      6402373705728000.0L,
      121645100408832000.0L,
      0.243290200817664e19L,
      0.5109094217170944e20L,
      0.112400072777760768e22L,
      0.2585201673888497664e23L,
      0.62044840173323943936e24L,
      0.15511210043330985984e26L,
      0.403291461126605635584e27L,
      0.10888869450418352160768e29L,
      0.304888344611713860501504e30L,
      0.8841761993739701954543616e31L,
      0.26525285981219105863630848e33L,
      0.822283865417792281772556288e34L,
      0.26313083693369353016721801216e36L,
      0.868331761881188649551819440128e37L,
      0.29523279903960414084761860964352e39L,
      0.103331479663861449296666513375232e41L,
      0.3719933267899012174679994481508352e42L,
      0.137637530912263450463159795815809024e44L,
      0.5230226174666011117600072241000742912e45L,
      0.203978820811974433586402817399028973568e47L,
      0.815915283247897734345611269596115894272e48L,
      0.3345252661316380710817006205344075166515e50L,
      0.1405006117752879898543142606244511569936e52L,
      0.6041526306337383563735513206851399750726e53L,
      0.265827157478844876804362581101461589032e55L,
      0.1196222208654801945619631614956577150644e57L,
      0.5502622159812088949850305428800254892962e58L,
      0.2586232415111681806429643551536119799692e60L,
      0.1241391559253607267086228904737337503852e62L,
      0.6082818640342675608722521633212953768876e63L,
      0.3041409320171337804361260816606476884438e65L,
      0.1551118753287382280224243016469303211063e67L,
      0.8065817517094387857166063685640376697529e68L,
      0.427488328406002556429801375338939964969e70L,
      0.2308436973392413804720927426830275810833e72L,
      0.1269640335365827592596510084756651695958e74L,
      0.7109985878048634518540456474637249497365e75L,
      0.4052691950487721675568060190543232213498e77L,
      0.2350561331282878571829474910515074683829e79L,
      0.1386831185456898357379390197203894063459e81L,
      0.8320987112741390144276341183223364380754e82L,
      0.507580213877224798800856812176625227226e84L,
      0.3146997326038793752565312235495076408801e86L,
      0.1982608315404440064116146708361898137545e88L,
      0.1268869321858841641034333893351614808029e90L,
      0.8247650592082470666723170306785496252186e91L,
      0.5443449390774430640037292402478427526443e93L,
      0.3647111091818868528824985909660546442717e95L,
      0.2480035542436830599600990418569171581047e97L,
      0.1711224524281413113724683388812728390923e99L,
      0.1197857166996989179607278372168909873646e101L,
      0.8504785885678623175211676442399260102886e102L,
      0.6123445837688608686152407038527467274078e104L,
      0.4470115461512684340891257138125051110077e106L,
      0.3307885441519386412259530282212537821457e108L,
      0.2480914081139539809194647711659403366093e110L,
      0.188549470166605025498793226086114655823e112L,
      0.1451830920282858696340707840863082849837e114L,
      0.1132428117820629783145752115873204622873e116L,
      0.8946182130782975286851441715398316520698e117L,
      0.7156945704626380229481153372318653216558e119L,
      0.5797126020747367985879734231578109105412e121L,
      0.4753643337012841748421382069894049466438e123L,
      0.3945523969720658651189747118012061057144e125L,
      0.3314240134565353266999387579130131288001e127L,
      0.2817104114380550276949479442260611594801e129L,
      0.2422709538367273238176552320344125971528e131L,
      0.210775729837952771721360051869938959523e133L,
      0.1854826422573984391147968456455462843802e135L,
      0.1650795516090846108121691926245361930984e137L,
      0.1485715964481761497309522733620825737886e139L,
      0.1352001527678402962551665687594951421476e141L,
      0.1243841405464130725547532432587355307758e143L,
      0.1156772507081641574759205162306240436215e145L,
      0.1087366156656743080273652852567866010042e147L,
      0.103299784882390592625997020993947270954e149L,
      0.9916779348709496892095714015418938011582e150L,
      0.9619275968248211985332842594956369871234e152L,
      0.942689044888324774562618574305724247381e154L,
      0.9332621544394415268169923885626670049072e156L,
      0.9332621544394415268169923885626670049072e158L,
      0.9425947759838359420851623124482936749562e160L,
      0.9614466715035126609268655586972595484554e162L,
      0.990290071648618040754671525458177334909e164L,
      0.1029901674514562762384858386476504428305e167L,
      0.1081396758240290900504101305800329649721e169L,
      0.1146280563734708354534347384148349428704e171L,
      0.1226520203196137939351751701038733888713e173L,
      0.132464181945182897449989183712183259981e175L,
      0.1443859583202493582204882102462797533793e177L,
      0.1588245541522742940425370312709077287172e179L,
      0.1762952551090244663872161047107075788761e181L,
      0.1974506857221074023536820372759924883413e183L,
      0.2231192748659813646596607021218715118256e185L,
      0.2543559733472187557120132004189335234812e187L,
      0.2925093693493015690688151804817735520034e189L,
      0.339310868445189820119825609358857320324e191L,
      0.396993716080872089540195962949863064779e193L,
      0.4684525849754290656574312362808384164393e195L,
      0.5574585761207605881323431711741977155627e197L,
      0.6689502913449127057588118054090372586753e199L,
      0.8094298525273443739681622845449350829971e201L,
      0.9875044200833601362411579871448208012564e203L,
      0.1214630436702532967576624324188129585545e206L,
      0.1506141741511140879795014161993280686076e208L,
      0.1882677176888926099743767702491600857595e210L,
      0.237217324288004688567714730513941708057e212L,
      0.3012660018457659544809977077527059692324e214L,
      0.3856204823625804217356770659234636406175e216L,
      0.4974504222477287440390234150412680963966e218L,
      0.6466855489220473672507304395536485253155e220L,
      0.8471580690878820510984568758152795681634e222L,
      0.1118248651196004307449963076076169029976e225L,
      0.1487270706090685728908450891181304809868e227L,
      0.1992942746161518876737324194182948445223e229L,
      0.269047270731805048359538766214698040105e231L,
      0.3659042881952548657689727220519893345429e233L,
      0.5012888748274991661034926292112253883237e235L,
      0.6917786472619488492228198283114910358867e237L,
      0.9615723196941089004197195613529725398826e239L,
      0.1346201247571752460587607385894161555836e242L,
      0.1898143759076170969428526414110767793728e244L,
      0.2695364137888162776588507508037290267094e246L,
      0.3854370717180072770521565736493325081944e248L,
      0.5550293832739304789551054660550388118e250L,
      0.80479260574719919448490292577980627711e252L,
      0.1174997204390910823947958271638517164581e255L,
      0.1727245890454638911203498659308620231933e257L,
      0.2556323917872865588581178015776757943262e259L,
      0.380892263763056972698595524350736933546e261L,
      0.571338395644585459047893286526105400319e263L,
      0.8627209774233240431623188626544191544816e265L,
      0.1311335885683452545606724671234717114812e268L,
      0.2006343905095682394778288746989117185662e270L,
      0.308976961384735088795856467036324046592e272L,
      0.4789142901463393876335775239063022722176e274L,
      0.7471062926282894447083809372938315446595e276L,
      0.1172956879426414428192158071551315525115e279L,
      0.1853271869493734796543609753051078529682e281L,
      0.2946702272495038326504339507351214862195e283L,
      0.4714723635992061322406943211761943779512e285L,
      0.7590705053947218729075178570936729485014e287L,
      0.1229694218739449434110178928491750176572e290L,
      0.2004401576545302577599591653441552787813e292L,
      0.3287218585534296227263330311644146572013e294L,
      0.5423910666131588774984495014212841843822e296L,
      0.9003691705778437366474261723593317460744e298L,
      0.1503616514864999040201201707840084015944e301L,
      0.2526075744973198387538018869171341146786e303L,
      0.4269068009004705274939251888899566538069e305L,
      0.7257415615307998967396728211129263114717e307L,
   }};

   return factorials[i];
}

template <>
struct max_factorial<long double>
{
   BOOST_STATIC_CONSTANT(unsigned, value = 170);
};

#ifdef BOOST_MATH_USE_FLOAT128

template <>
inline BOOST_MATH_FLOAT128_TYPE unchecked_factorial<BOOST_MATH_FLOAT128_TYPE>(unsigned i)
{
   static const boost::array<BOOST_MATH_FLOAT128_TYPE, 171> factorials = { {
      1,
      1,
      2,
      6,
      24,
      120,
      720,
      5040,
      40320,
      362880.0Q,
      3628800.0Q,
      39916800.0Q,
      479001600.0Q,
      6227020800.0Q,
      87178291200.0Q,
      1307674368000.0Q,
      20922789888000.0Q,
      355687428096000.0Q,
      6402373705728000.0Q,
      121645100408832000.0Q,
      0.243290200817664e19Q,
      0.5109094217170944e20Q,
      0.112400072777760768e22Q,
      0.2585201673888497664e23Q,
      0.62044840173323943936e24Q,
      0.15511210043330985984e26Q,
      0.403291461126605635584e27Q,
      0.10888869450418352160768e29Q,
      0.304888344611713860501504e30Q,
      0.8841761993739701954543616e31Q,
      0.26525285981219105863630848e33Q,
      0.822283865417792281772556288e34Q,
      0.26313083693369353016721801216e36Q,
      0.868331761881188649551819440128e37Q,
      0.29523279903960414084761860964352e39Q,
      0.103331479663861449296666513375232e41Q,
      0.3719933267899012174679994481508352e42Q,
      0.137637530912263450463159795815809024e44Q,
      0.5230226174666011117600072241000742912e45Q,
      0.203978820811974433586402817399028973568e47Q,
      0.815915283247897734345611269596115894272e48Q,
      0.3345252661316380710817006205344075166515e50Q,
      0.1405006117752879898543142606244511569936e52Q,
      0.6041526306337383563735513206851399750726e53Q,
      0.265827157478844876804362581101461589032e55Q,
      0.1196222208654801945619631614956577150644e57Q,
      0.5502622159812088949850305428800254892962e58Q,
      0.2586232415111681806429643551536119799692e60Q,
      0.1241391559253607267086228904737337503852e62Q,
      0.6082818640342675608722521633212953768876e63Q,
      0.3041409320171337804361260816606476884438e65Q,
      0.1551118753287382280224243016469303211063e67Q,
      0.8065817517094387857166063685640376697529e68Q,
      0.427488328406002556429801375338939964969e70Q,
      0.2308436973392413804720927426830275810833e72Q,
      0.1269640335365827592596510084756651695958e74Q,
      0.7109985878048634518540456474637249497365e75Q,
      0.4052691950487721675568060190543232213498e77Q,
      0.2350561331282878571829474910515074683829e79Q,
      0.1386831185456898357379390197203894063459e81Q,
      0.8320987112741390144276341183223364380754e82Q,
      0.507580213877224798800856812176625227226e84Q,
      0.3146997326038793752565312235495076408801e86Q,
      0.1982608315404440064116146708361898137545e88Q,
      0.1268869321858841641034333893351614808029e90Q,
      0.8247650592082470666723170306785496252186e91Q,
      0.5443449390774430640037292402478427526443e93Q,
      0.3647111091818868528824985909660546442717e95Q,
      0.2480035542436830599600990418569171581047e97Q,
      0.1711224524281413113724683388812728390923e99Q,
      0.1197857166996989179607278372168909873646e101Q,
      0.8504785885678623175211676442399260102886e102Q,
      0.6123445837688608686152407038527467274078e104Q,
      0.4470115461512684340891257138125051110077e106Q,
      0.3307885441519386412259530282212537821457e108Q,
      0.2480914081139539809194647711659403366093e110Q,
      0.188549470166605025498793226086114655823e112Q,
      0.1451830920282858696340707840863082849837e114Q,
      0.1132428117820629783145752115873204622873e116Q,
      0.8946182130782975286851441715398316520698e117Q,
      0.7156945704626380229481153372318653216558e119Q,
      0.5797126020747367985879734231578109105412e121Q,
      0.4753643337012841748421382069894049466438e123Q,
      0.3945523969720658651189747118012061057144e125Q,
      0.3314240134565353266999387579130131288001e127Q,
      0.2817104114380550276949479442260611594801e129Q,
      0.2422709538367273238176552320344125971528e131Q,
      0.210775729837952771721360051869938959523e133Q,
      0.1854826422573984391147968456455462843802e135Q,
      0.1650795516090846108121691926245361930984e137Q,
      0.1485715964481761497309522733620825737886e139Q,
      0.1352001527678402962551665687594951421476e141Q,
      0.1243841405464130725547532432587355307758e143Q,
      0.1156772507081641574759205162306240436215e145Q,
      0.1087366156656743080273652852567866010042e147Q,
      0.103299784882390592625997020993947270954e149Q,
      0.9916779348709496892095714015418938011582e150Q,
      0.9619275968248211985332842594956369871234e152Q,
      0.942689044888324774562618574305724247381e154Q,
      0.9332621544394415268169923885626670049072e156Q,
      0.9332621544394415268169923885626670049072e158Q,
      0.9425947759838359420851623124482936749562e160Q,
      0.9614466715035126609268655586972595484554e162Q,
      0.990290071648618040754671525458177334909e164Q,
      0.1029901674514562762384858386476504428305e167Q,
      0.1081396758240290900504101305800329649721e169Q,
      0.1146280563734708354534347384148349428704e171Q,
      0.1226520203196137939351751701038733888713e173Q,
      0.132464181945182897449989183712183259981e175Q,
      0.1443859583202493582204882102462797533793e177Q,
      0.1588245541522742940425370312709077287172e179Q,
      0.1762952551090244663872161047107075788761e181Q,
      0.1974506857221074023536820372759924883413e183Q,
      0.2231192748659813646596607021218715118256e185Q,
      0.2543559733472187557120132004189335234812e187Q,
      0.2925093693493015690688151804817735520034e189Q,
      0.339310868445189820119825609358857320324e191Q,
      0.396993716080872089540195962949863064779e193Q,
      0.4684525849754290656574312362808384164393e195Q,
      0.5574585761207605881323431711741977155627e197Q,
      0.6689502913449127057588118054090372586753e199Q,
      0.8094298525273443739681622845449350829971e201Q,
      0.9875044200833601362411579871448208012564e203Q,
      0.1214630436702532967576624324188129585545e206Q,
      0.1506141741511140879795014161993280686076e208Q,
      0.1882677176888926099743767702491600857595e210Q,
      0.237217324288004688567714730513941708057e212Q,
      0.3012660018457659544809977077527059692324e214Q,
      0.3856204823625804217356770659234636406175e216Q,
      0.4974504222477287440390234150412680963966e218Q,
      0.6466855489220473672507304395536485253155e220Q,
      0.8471580690878820510984568758152795681634e222Q,
      0.1118248651196004307449963076076169029976e225Q,
      0.1487270706090685728908450891181304809868e227Q,
      0.1992942746161518876737324194182948445223e229Q,
      0.269047270731805048359538766214698040105e231Q,
      0.3659042881952548657689727220519893345429e233Q,
      0.5012888748274991661034926292112253883237e235Q,
      0.6917786472619488492228198283114910358867e237Q,
      0.9615723196941089004197195613529725398826e239Q,
      0.1346201247571752460587607385894161555836e242Q,
      0.1898143759076170969428526414110767793728e244Q,
      0.2695364137888162776588507508037290267094e246Q,
      0.3854370717180072770521565736493325081944e248Q,
      0.5550293832739304789551054660550388118e250Q,
      0.80479260574719919448490292577980627711e252Q,
      0.1174997204390910823947958271638517164581e255Q,
      0.1727245890454638911203498659308620231933e257Q,
      0.2556323917872865588581178015776757943262e259Q,
      0.380892263763056972698595524350736933546e261Q,
      0.571338395644585459047893286526105400319e263Q,
      0.8627209774233240431623188626544191544816e265Q,
      0.1311335885683452545606724671234717114812e268Q,
      0.2006343905095682394778288746989117185662e270Q,
      0.308976961384735088795856467036324046592e272Q,
      0.4789142901463393876335775239063022722176e274Q,
      0.7471062926282894447083809372938315446595e276Q,
      0.1172956879426414428192158071551315525115e279Q,
      0.1853271869493734796543609753051078529682e281Q,
      0.2946702272495038326504339507351214862195e283Q,
      0.4714723635992061322406943211761943779512e285Q,
      0.7590705053947218729075178570936729485014e287Q,
      0.1229694218739449434110178928491750176572e290Q,
      0.2004401576545302577599591653441552787813e292Q,
      0.3287218585534296227263330311644146572013e294Q,
      0.5423910666131588774984495014212841843822e296Q,
      0.9003691705778437366474261723593317460744e298Q,
      0.1503616514864999040201201707840084015944e301Q,
      0.2526075744973198387538018869171341146786e303Q,
      0.4269068009004705274939251888899566538069e305Q,
      0.7257415615307998967396728211129263114717e307Q,
   } };

   return factorials[i];
}

template <>
struct max_factorial<BOOST_MATH_FLOAT128_TYPE>
{
   BOOST_STATIC_CONSTANT(unsigned, value = 170);
};

#endif

template <>
inline double unchecked_factorial<double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(double))
{
   return static_cast<double>(boost::math::unchecked_factorial<long double>(i));
}

template <>
struct max_factorial<double>
{
   BOOST_STATIC_CONSTANT(unsigned,
      value = ::boost::math::max_factorial<long double>::value);
};

#ifndef BOOST_MATH_NO_LEXICAL_CAST

template <class T>
struct unchecked_factorial_initializer
{
   struct init
   {
      init()
      {
         boost::math::unchecked_factorial<T>(3);
      }
      void force_instantiate()const {}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

template <class T>
const typename unchecked_factorial_initializer<T>::init unchecked_factorial_initializer<T>::initializer;

template <class T, int N>
inline T unchecked_factorial_imp(unsigned i, const mpl::int_<N>&)
{
   BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
   // factorial<unsigned int>(n) is not implemented
   // because it would overflow integral type T for too small n
   // to be useful. Use instead a floating-point type,
   // and convert to an unsigned type if essential, for example:
   // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n));
   // See factorial documentation for more detail.

   unchecked_factorial_initializer<T>::force_instantiate();

   static const boost::array<T, 101> factorials = {{
      T(boost::math::tools::convert_from_string<T>("1")),
      T(boost::math::tools::convert_from_string<T>("1")),
      T(boost::math::tools::convert_from_string<T>("2")),
      T(boost::math::tools::convert_from_string<T>("6")),
      T(boost::math::tools::convert_from_string<T>("24")),
      T(boost::math::tools::convert_from_string<T>("120")),
      T(boost::math::tools::convert_from_string<T>("720")),
      T(boost::math::tools::convert_from_string<T>("5040")),
      T(boost::math::tools::convert_from_string<T>("40320")),
      T(boost::math::tools::convert_from_string<T>("362880")),
      T(boost::math::tools::convert_from_string<T>("3628800")),
      T(boost::math::tools::convert_from_string<T>("39916800")),
      T(boost::math::tools::convert_from_string<T>("479001600")),
      T(boost::math::tools::convert_from_string<T>("6227020800")),
      T(boost::math::tools::convert_from_string<T>("87178291200")),
      T(boost::math::tools::convert_from_string<T>("1307674368000")),
      T(boost::math::tools::convert_from_string<T>("20922789888000")),
      T(boost::math::tools::convert_from_string<T>("355687428096000")),
      T(boost::math::tools::convert_from_string<T>("6402373705728000")),
      T(boost::math::tools::convert_from_string<T>("121645100408832000")),
      T(boost::math::tools::convert_from_string<T>("2432902008176640000")),
      T(boost::math::tools::convert_from_string<T>("51090942171709440000")),
      T(boost::math::tools::convert_from_string<T>("1124000727777607680000")),
      T(boost::math::tools::convert_from_string<T>("25852016738884976640000")),
      T(boost::math::tools::convert_from_string<T>("620448401733239439360000")),
      T(boost::math::tools::convert_from_string<T>("15511210043330985984000000")),
      T(boost::math::tools::convert_from_string<T>("403291461126605635584000000")),
      T(boost::math::tools::convert_from_string<T>("10888869450418352160768000000")),
      T(boost::math::tools::convert_from_string<T>("304888344611713860501504000000")),
      T(boost::math::tools::convert_from_string<T>("8841761993739701954543616000000")),
      T(boost::math::tools::convert_from_string<T>("265252859812191058636308480000000")),
      T(boost::math::tools::convert_from_string<T>("8222838654177922817725562880000000")),
      T(boost::math::tools::convert_from_string<T>("263130836933693530167218012160000000")),
      T(boost::math::tools::convert_from_string<T>("8683317618811886495518194401280000000")),
      T(boost::math::tools::convert_from_string<T>("295232799039604140847618609643520000000")),
      T(boost::math::tools::convert_from_string<T>("10333147966386144929666651337523200000000")),
      T(boost::math::tools::convert_from_string<T>("371993326789901217467999448150835200000000")),
      T(boost::math::tools::convert_from_string<T>("13763753091226345046315979581580902400000000")),
      T(boost::math::tools::convert_from_string<T>("523022617466601111760007224100074291200000000")),
      T(boost::math::tools::convert_from_string<T>("20397882081197443358640281739902897356800000000")),
      T(boost::math::tools::convert_from_string<T>("815915283247897734345611269596115894272000000000")),
      T(boost::math::tools::convert_from_string<T>("33452526613163807108170062053440751665152000000000")),
      T(boost::math::tools::convert_from_string<T>("1405006117752879898543142606244511569936384000000000")),
      T(boost::math::tools::convert_from_string<T>("60415263063373835637355132068513997507264512000000000")),
      T(boost::math::tools::convert_from_string<T>("2658271574788448768043625811014615890319638528000000000")),
      T(boost::math::tools::convert_from_string<T>("119622220865480194561963161495657715064383733760000000000")),
      T(boost::math::tools::convert_from_string<T>("5502622159812088949850305428800254892961651752960000000000")),
      T(boost::math::tools::convert_from_string<T>("258623241511168180642964355153611979969197632389120000000000")),
      T(boost::math::tools::convert_from_string<T>("12413915592536072670862289047373375038521486354677760000000000")),
      T(boost::math::tools::convert_from_string<T>("608281864034267560872252163321295376887552831379210240000000000")),
      T(boost::math::tools::convert_from_string<T>("30414093201713378043612608166064768844377641568960512000000000000")),
      T(boost::math::tools::convert_from_string<T>("1551118753287382280224243016469303211063259720016986112000000000000")),
      T(boost::math::tools::convert_from_string<T>("80658175170943878571660636856403766975289505440883277824000000000000")),
      T(boost::math::tools::convert_from_string<T>("4274883284060025564298013753389399649690343788366813724672000000000000")),
      T(boost::math::tools::convert_from_string<T>("230843697339241380472092742683027581083278564571807941132288000000000000")),
      T(boost::math::tools::convert_from_string<T>("12696403353658275925965100847566516959580321051449436762275840000000000000")),
      T(boost::math::tools::convert_from_string<T>("710998587804863451854045647463724949736497978881168458687447040000000000000")),
      T(boost::math::tools::convert_from_string<T>("40526919504877216755680601905432322134980384796226602145184481280000000000000")),
      T(boost::math::tools::convert_from_string<T>("2350561331282878571829474910515074683828862318181142924420699914240000000000000")),
      T(boost::math::tools::convert_from_string<T>("138683118545689835737939019720389406345902876772687432540821294940160000000000000")),
      T(boost::math::tools::convert_from_string<T>("8320987112741390144276341183223364380754172606361245952449277696409600000000000000")),
      T(boost::math::tools::convert_from_string<T>("507580213877224798800856812176625227226004528988036003099405939480985600000000000000")),
      T(boost::math::tools::convert_from_string<T>("31469973260387937525653122354950764088012280797258232192163168247821107200000000000000")),
      T(boost::math::tools::convert_from_string<T>("1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000")),
      T(boost::math::tools::convert_from_string<T>("126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000")),
      T(boost::math::tools::convert_from_string<T>("8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000")),
      T(boost::math::tools::convert_from_string<T>("544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000")),
      T(boost::math::tools::convert_from_string<T>("36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000")),
      T(boost::math::tools::convert_from_string<T>("2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000")),
      T(boost::math::tools::convert_from_string<T>("171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000")),
      T(boost::math::tools::convert_from_string<T>("11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000")),
      T(boost::math::tools::convert_from_string<T>("850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000")),
      T(boost::math::tools::convert_from_string<T>("61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000")),
      T(boost::math::tools::convert_from_string<T>("4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000")),
      T(boost::math::tools::convert_from_string<T>("330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000")),
      T(boost::math::tools::convert_from_string<T>("24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000")),
      T(boost::math::tools::convert_from_string<T>("93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000")),
   }};

   return factorials[i];
}

template <class T>
inline T unchecked_factorial_imp(unsigned i, const mpl::int_<0>&)
{
   BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
   // factorial<unsigned int>(n) is not implemented
   // because it would overflow integral type T for too small n
   // to be useful. Use instead a floating-point type,
   // and convert to an unsigned type if essential, for example:
   // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n));
   // See factorial documentation for more detail.
#ifdef BOOST_NO_CXX11_THREAD_LOCAL
   unchecked_factorial_initializer<T>::force_instantiate();
#endif
   static const char* const factorial_strings[] = {
         "1",
         "1",
         "2",
         "6",
         "24",
         "120",
         "720",
         "5040",
         "40320",
         "362880",
         "3628800",
         "39916800",
         "479001600",
         "6227020800",
         "87178291200",
         "1307674368000",
         "20922789888000",
         "355687428096000",
         "6402373705728000",
         "121645100408832000",
         "2432902008176640000",
         "51090942171709440000",
         "1124000727777607680000",
         "25852016738884976640000",
         "620448401733239439360000",
         "15511210043330985984000000",
         "403291461126605635584000000",
         "10888869450418352160768000000",
         "304888344611713860501504000000",
         "8841761993739701954543616000000",
         "265252859812191058636308480000000",
         "8222838654177922817725562880000000",
         "263130836933693530167218012160000000",
         "8683317618811886495518194401280000000",
         "295232799039604140847618609643520000000",
         "10333147966386144929666651337523200000000",
         "371993326789901217467999448150835200000000",
         "13763753091226345046315979581580902400000000",
         "523022617466601111760007224100074291200000000",
         "20397882081197443358640281739902897356800000000",
         "815915283247897734345611269596115894272000000000",
         "33452526613163807108170062053440751665152000000000",
         "1405006117752879898543142606244511569936384000000000",
         "60415263063373835637355132068513997507264512000000000",
         "2658271574788448768043625811014615890319638528000000000",
         "119622220865480194561963161495657715064383733760000000000",
         "5502622159812088949850305428800254892961651752960000000000",
         "258623241511168180642964355153611979969197632389120000000000",
         "12413915592536072670862289047373375038521486354677760000000000",
         "608281864034267560872252163321295376887552831379210240000000000",
         "30414093201713378043612608166064768844377641568960512000000000000",
         "1551118753287382280224243016469303211063259720016986112000000000000",
         "80658175170943878571660636856403766975289505440883277824000000000000",
         "4274883284060025564298013753389399649690343788366813724672000000000000",
         "230843697339241380472092742683027581083278564571807941132288000000000000",
         "12696403353658275925965100847566516959580321051449436762275840000000000000",
         "710998587804863451854045647463724949736497978881168458687447040000000000000",
         "40526919504877216755680601905432322134980384796226602145184481280000000000000",
         "2350561331282878571829474910515074683828862318181142924420699914240000000000000",
         "138683118545689835737939019720389406345902876772687432540821294940160000000000000",
         "8320987112741390144276341183223364380754172606361245952449277696409600000000000000",
         "507580213877224798800856812176625227226004528988036003099405939480985600000000000000",
         "31469973260387937525653122354950764088012280797258232192163168247821107200000000000000",
         "1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000",
         "126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000",
         "8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000",
         "544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000",
         "36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000",
         "2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000",
         "171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000",
         "11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000",
         "850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000",
         "61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000",
         "4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000",
         "330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000",
         "24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000",
         "1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000",
         "145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000",
         "11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000",
         "894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000",
         "71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000",
         "5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000",
         "475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000",
         "39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000",
         "3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000",
         "281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000",
         "24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000",
         "2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000",
         "185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000",
         "16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000",
         "1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000",
         "135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000",
         "12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000",
         "1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000",
         "108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000",
         "10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000",
         "991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000",
         "96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000",
         "9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000",
         "933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000",
         "93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000",
      };

      static BOOST_MATH_THREAD_LOCAL T factorials[sizeof(factorial_strings) / sizeof(factorial_strings[0])];
      static BOOST_MATH_THREAD_LOCAL int digits = 0;

      int current_digits = boost::math::tools::digits<T>();

      if(digits != current_digits)
      {
         digits = current_digits;
         for(unsigned k = 0; k < sizeof(factorials) / sizeof(factorials[0]); ++k)
            factorials[k] = static_cast<T>(boost::math::tools::convert_from_string<T>(factorial_strings[k]));
      }

   return factorials[i];
}

template <class T>
inline T unchecked_factorial(unsigned i)
{
   typedef typename boost::math::policies::precision<T, boost::math::policies::policy<> >::type tag_type;
   return unchecked_factorial_imp<T>(i, tag_type());
}

template <class T>
struct max_factorial
{
   BOOST_STATIC_CONSTANT(unsigned, value = 100);
};

#else // BOOST_MATH_NO_LEXICAL_CAST

template <class T>
inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
{
   return 1;
}

template <class T>
struct max_factorial
{
   BOOST_STATIC_CONSTANT(unsigned, value = 0);
};

#endif

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
template <class T>
const unsigned max_factorial<T>::value;
#endif

} // namespace math
} // namespace boost

#endif // BOOST_MATH_SP_UC_FACTORIALS_HPP

license: mit


code:
//  (C) Copyright John Maddock 2006.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_SF_DIGAMMA_HPP
#define BOOST_MATH_SF_DIGAMMA_HPP

#ifdef _MSC_VER
#pragma once
#pragma warning(push)
#pragma warning(disable:4702) // Unreachable code (release mode only warning)
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/mpl/comparison.hpp>
#include <boost/math/tools/big_constant.hpp>

namespace boost{
namespace math{
namespace detail{
//
// Begin by defining the smallest value for which it is safe to
// use the asymptotic expansion for digamma:
//
inline unsigned digamma_large_lim(const mpl::int_<0>*)
{  return 20;  }
inline unsigned digamma_large_lim(const mpl::int_<113>*)
{  return 20;  }
inline unsigned digamma_large_lim(const void*)
{  return 10;  }
//
// Implementations of the asymptotic expansion come next,
// the coefficients of the series have been evaluated
// in advance at high precision, and the series truncated
// at the first term that's too small to effect the result.
// Note that the series becomes divergent after a while
// so truncation is very important.
//
// This first one gives 34-digit precision for x >= 20:
//
template <class T>
inline T digamma_imp_large(T x, const mpl::int_<113>*)
{
   BOOST_MATH_STD_USING // ADL of std functions.
   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.0083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.003968253968253968253968253968253968253968253968254),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.0041666666666666666666666666666666666666666666666667),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.0075757575757575757575757575757575757575757575757576),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.021092796092796092796092796092796092796092796092796),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.44325980392156862745098039215686274509803921568627),
      BOOST_MATH_BIG_CONSTANT(T, 113, 3.0539543302701197438039543302701197438039543302701),
      BOOST_MATH_BIG_CONSTANT(T, 113, -26.456212121212121212121212121212121212121212121212),
      BOOST_MATH_BIG_CONSTANT(T, 113, 281.4601449275362318840579710144927536231884057971),
      BOOST_MATH_BIG_CONSTANT(T, 113, -3607.510546398046398046398046398046398046398046398),
      BOOST_MATH_BIG_CONSTANT(T, 113, 54827.583333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 113, -974936.82385057471264367816091954022988505747126437),
      BOOST_MATH_BIG_CONSTANT(T, 113, 20052695.796688078946143462272494530559046688078946),
      BOOST_MATH_BIG_CONSTANT(T, 113, -472384867.72162990196078431372549019607843137254902),
      BOOST_MATH_BIG_CONSTANT(T, 113, 12635724795.916666666666666666666666666666666666667)
   };
   x -= 1;
   T result = log(x);
   result += 1 / (2 * x);
   T z = 1 / (x*x);
   result -= z * tools::evaluate_polynomial(P, z);
   return result;
}
//
// 19-digit precision for x >= 10:
//
template <class T>
inline T digamma_imp_large(T x, const mpl::int_<64>*)
{
   BOOST_MATH_STD_USING // ADL of std functions.
   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.0083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.003968253968253968253968253968253968253968253968254),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.0041666666666666666666666666666666666666666666666667),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.0075757575757575757575757575757575757575757575757576),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.021092796092796092796092796092796092796092796092796),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.44325980392156862745098039215686274509803921568627),
      BOOST_MATH_BIG_CONSTANT(T, 64, 3.0539543302701197438039543302701197438039543302701),
      BOOST_MATH_BIG_CONSTANT(T, 64, -26.456212121212121212121212121212121212121212121212),
      BOOST_MATH_BIG_CONSTANT(T, 64, 281.4601449275362318840579710144927536231884057971),
   };
   x -= 1;
   T result = log(x);
   result += 1 / (2 * x);
   T z = 1 / (x*x);
   result -= z * tools::evaluate_polynomial(P, z);
   return result;
}
//
// 17-digit precision for x >= 10:
//
template <class T>
inline T digamma_imp_large(T x, const mpl::int_<53>*)
{
   BOOST_MATH_STD_USING // ADL of std functions.
   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.0083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.003968253968253968253968253968253968253968253968254),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.0041666666666666666666666666666666666666666666666667),
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.0075757575757575757575757575757575757575757575757576),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.021092796092796092796092796092796092796092796092796),
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.44325980392156862745098039215686274509803921568627)
   };
   x -= 1;
   T result = log(x);
   result += 1 / (2 * x);
   T z = 1 / (x*x);
   result -= z * tools::evaluate_polynomial(P, z);
   return result;
}
//
// 9-digit precision for x >= 10:
//
template <class T>
inline T digamma_imp_large(T x, const mpl::int_<24>*)
{
   BOOST_MATH_STD_USING // ADL of std functions.
   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 24, 0.083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 24, -0.0083333333333333333333333333333333333333333333333333),
      BOOST_MATH_BIG_CONSTANT(T, 24, 0.003968253968253968253968253968253968253968253968254)
   };
   x -= 1;
   T result = log(x);
   result += 1 / (2 * x);
   T z = 1 / (x*x);
   result -= z * tools::evaluate_polynomial(P, z);
   return result;
}
//
// Fully generic asymptotic expansion in terms of Bernoulli numbers, see:
// http://functions.wolfram.com/06.14.06.0012.01
//
template <class T>
struct digamma_series_func
{
private:
   int k;
   T xx;
   T term;
public:
   digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {}
   T operator()()
   {
      T result = term * boost::math::bernoulli_b2n<T>(k) / (2 * k);
      term /= xx;
      ++k;
      return result;
   }
   typedef T result_type;
};

template <class T, class Policy>
inline T digamma_imp_large(T x, const Policy& pol, const mpl::int_<0>*)
{
   BOOST_MATH_STD_USING
   digamma_series_func<T> s(x);
   T result = log(x) - 1 / (2 * x);
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
   result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, -result);
   result = -result;
   policies::check_series_iterations<T>("boost::math::digamma<%1%>(%1%)", max_iter, pol);
   return result;
}
//
// Now follow rational approximations over the range [1,2].
//
// 35-digit precision:
//
template <class T>
T digamma_imp_1_2(T x, const mpl::int_<113>*)
{
   //
   // Now the approximation, we use the form:
   //
   // digamma(x) = (x - root) * (Y + R(x-1))
   //
   // Where root is the location of the positive root of digamma,
   // Y is a constant, and R is optimised for low absolute error
   // compared to Y.
   //
   // Max error found at 128-bit long double precision:  5.541e-35
   // Maximum Deviation Found (approximation error):     1.965e-35
   //
   static const float Y = 0.99558162689208984375F;

   static const T root1 = T(1569415565) / 1073741824uL;
   static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
   static const T root3 = ((T(111616537) / 1073741824uL) / 1073741824uL) / 1073741824uL;
   static const T root4 = (((T(503992070) / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL;
   static const T root5 = BOOST_MATH_BIG_CONSTANT(T, 113, 0.52112228569249997894452490385577338504019838794544e-36);

   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.25479851061131551526977464225335883769),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.18684290534374944114622235683619897417),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.80360876047931768958995775910991929922),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.67227342794829064330498117008564270136),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.26569010991230617151285010695543858005),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.05775672694575986971640757748003553385),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.0071432147823164975485922555833274240665),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.00048740753910766168912364555706064993274),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.16454996865214115723416538844975174761e-4),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.20327832297631728077731148515093164955e-6)
   };
   static const T Q[] = {
      BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
      BOOST_MATH_BIG_CONSTANT(T, 113, 2.6210924610812025425088411043163287646),
      BOOST_MATH_BIG_CONSTANT(T, 113, 2.6850757078559596612621337395886392594),
      BOOST_MATH_BIG_CONSTANT(T, 113, 1.4320913706209965531250495490639289418),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.4410872083455009362557012239501953402),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.081385727399251729505165509278152487225),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089478633066857163432104815183858149496),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.00055861622855066424871506755481997374154),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.1760168552357342401304462967950178554e-4),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.20585454493572473724556649516040874384e-6),
      BOOST_MATH_BIG_CONSTANT(T, 113, -0.90745971844439990284514121823069162795e-11),
      BOOST_MATH_BIG_CONSTANT(T, 113, 0.48857673606545846774761343500033283272e-13),
   };
   T g = x - root1;
   g -= root2;
   g -= root3;
   g -= root4;
   g -= root5;
   T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
   T result = g * Y + g * r;

   return result;
}
//
// 19-digit precision:
//
template <class T>
T digamma_imp_1_2(T x, const mpl::int_<64>*)
{
   //
   // Now the approximation, we use the form:
   //
   // digamma(x) = (x - root) * (Y + R(x-1))
   //
   // Where root is the location of the positive root of digamma,
   // Y is a constant, and R is optimised for low absolute error
   // compared to Y.
   //
   // Max error found at 80-bit long double precision:   5.016e-20
   // Maximum Deviation Found (approximation error):     3.575e-20
   //
   static const float Y = 0.99558162689208984375F;

   static const T root1 = T(1569415565) / 1073741824uL;
   static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
   static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.9016312093258695918615325266959189453125e-19);

   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.254798510611315515235),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.314628554532916496608),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.665836341559876230295),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.314767657147375752913),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.0541156266153505273939),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.00289268368333918761452)
   };
   static const T Q[] = {
      BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
      BOOST_MATH_BIG_CONSTANT(T, 64, 2.1195759927055347547),
      BOOST_MATH_BIG_CONSTANT(T, 64, 1.54350554664961128724),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.486986018231042975162),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.0660481487173569812846),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.00298999662592323990972),
      BOOST_MATH_BIG_CONSTANT(T, 64, -0.165079794012604905639e-5),
      BOOST_MATH_BIG_CONSTANT(T, 64, 0.317940243105952177571e-7)
   };
   T g = x - root1;
   g -= root2;
   g -= root3;
   T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
   T result = g * Y + g * r;

   return result;
}
//
// 18-digit precision:
//
template <class T>
T digamma_imp_1_2(T x, const mpl::int_<53>*)
{
   //
   // Now the approximation, we use the form:
   //
   // digamma(x) = (x - root) * (Y + R(x-1))
   //
   // Where root is the location of the positive root of digamma,
   // Y is a constant, and R is optimised for low absolute error
   // compared to Y.
   //
   // Maximum Deviation Found:               1.466e-18
   // At double precision, max error found:  2.452e-17
   //
   static const float Y = 0.99558162689208984F;

   static const T root1 = T(1569415565) / 1073741824uL;
   static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
   static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.9016312093258695918615325266959189453125e-19);

   static const T P[] = {
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.25479851061131551),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.32555031186804491),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.65031853770896507),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.28919126444774784),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.045251321448739056),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.0020713321167745952)
   };
   static const T Q[] = {
      BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
      BOOST_MATH_BIG_CONSTANT(T, 53, 2.0767117023730469),
      BOOST_MATH_BIG_CONSTANT(T, 53, 1.4606242909763515),
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.43593529692665969),
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.054151797245674225),
      BOOST_MATH_BIG_CONSTANT(T, 53, 0.0021284987017821144),
      BOOST_MATH_BIG_CONSTANT(T, 53, -0.55789841321675513e-6)
   };
   T g = x - root1;
   g -= root2;
   g -= root3;
   T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
   T result = g * Y + g * r;

   return result;
}
//
// 9-digit precision:
//
template <class T>
inline T digamma_imp_1_2(T x, const mpl::int_<24>*)
{
   //
   // Now the approximation, we use the form:
   //
   // digamma(x) = (x - root) * (Y + R(x-1))
   //
   // Where root is the location of the positive root of digamma,
   // Y is a constant, and R is optimised for low absolute error
   // compared to Y.
   //
   // Maximum Deviation Found:              3.388e-010
   // At float precision, max error found:  2.008725e-008
   //
   static const float Y = 0.99558162689208984f;
   static const T root = 1532632.0f / 1048576;
   static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L);
   static const T P[] = {
      0.25479851023250261e0f,
      -0.44981331915268368e0f,
      -0.43916936919946835e0f,
      -0.61041765350579073e-1f
   };
   static const T Q[] = {
      0.1e1,
      0.15890202430554952e1f,
      0.65341249856146947e0f,
      0.63851690523355715e-1f
   };
   T g = x - root;
   g -= root_minor;
   T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
   T result = g * Y + g * r;

   return result;
}

template <class T, class Tag, class Policy>
T digamma_imp(T x, const Tag* t, const Policy& pol)
{
   //
   // This handles reflection of negative arguments, and all our
   // error handling, then forwards to the T-specific approximation.
   //
   BOOST_MATH_STD_USING // ADL of std functions.

   T result = 0;
   //
   // Check for negative arguments and use reflection:
   //
   if(x <= -1)
   {
      // Reflect:
      x = 1 - x;
      // Argument reduction for tan:
      T remainder = x - floor(x);
      // Shift to negative if > 0.5:
      if(remainder > 0.5)
      {
         remainder -= 1;
      }
      //
      // check for evaluation at a negative pole:
      //
      if(remainder == 0)
      {
         return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
      }
      result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
   }
   if(x == 0)
      return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
   //
   // If we're above the lower-limit for the
   // asymptotic expansion then use it:
   //
   if(x >= digamma_large_lim(t))
   {
      result += digamma_imp_large(x, t);
   }
   else
   {
      //
      // If x > 2 reduce to the interval [1,2]:
      //
      while(x > 2)
      {
         x -= 1;
         result += 1/x;
      }
      //
      // If x < 1 use recurrance to shift to > 1:
      //
      while(x < 1)
      {
         result -= 1/x;
         x += 1;
      }
      result += digamma_imp_1_2(x, t);
   }
   return result;
}

template <class T, class Policy>
T digamma_imp(T x, const mpl::int_<0>* t, const Policy& pol)
{
   //
   // This handles reflection of negative arguments, and all our
   // error handling, then forwards to the T-specific approximation.
   //
   BOOST_MATH_STD_USING // ADL of std functions.

   T result = 0;
   //
   // Check for negative arguments and use reflection:
   //
   if(x <= -1)
   {
      // Reflect:
      x = 1 - x;
      // Argument reduction for tan:
      T remainder = x - floor(x);
      // Shift to negative if > 0.5:
      if(remainder > 0.5)
      {
         remainder -= 1;
      }
      //
      // check for evaluation at a negative pole:
      //
      if(remainder == 0)
      {
         return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1 - x), pol);
      }
      result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
   }
   if(x == 0)
      return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
   //
   // If we're above the lower-limit for the
   // asymptotic expansion then use it, the
   // limit is a linear interpolation with
   // limit = 10 at 50 bit precision and
   // limit = 250 at 1000 bit precision.
   //
   int lim = 10 + ((tools::digits<T>() - 50) * 240L) / 950;
   T two_x = ldexp(x, 1);
   if(x >= lim)
   {
      result += digamma_imp_large(x, pol, t);
   }
   else if(floor(x) == x)
   {
      //
      // Special case for integer arguments, see
      // http://functions.wolfram.com/06.14.03.0001.01
      //
      result = -constants::euler<T, Policy>();
      T val = 1;
      while(val < x)
      {
         result += 1 / val;
         val += 1;
      }
   }
   else if(floor(two_x) == two_x)
   {
      //
      // Special case for half integer arguments, see:
      // http://functions.wolfram.com/06.14.03.0007.01
      //
      result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>();
      int n = itrunc(x);
      if(n)
      {
         for(int k = 1; k < n; ++k)
            result += 1 / T(k);
         for(int k = n; k <= 2 * n - 1; ++k)
            result += 2 / T(k);
      }
   }
   else
   {
      //
      // Rescale so we can use the asymptotic expansion:
      //
      while(x < lim)
      {
         result -= 1 / x;
         x += 1;
      }
      result += digamma_imp_large(x, pol, t);
   }
   return result;
}
//
// Initializer: ensure all our constants are initialized prior to the first call of main:
//
template <class T, class Policy>
struct digamma_initializer
{
   struct init
   {
      init()
      {
         typedef typename policies::precision<T, Policy>::type precision_type;
         do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>());
      }
      void do_init(const mpl::true_&)
      {
         boost::math::digamma(T(1.5), Policy());
         boost::math::digamma(T(500), Policy());
      }
      void do_init(const mpl::false_&){}
      void force_instantiate()const{}
   };
   static const init initializer;
   static void force_instantiate()
   {
      initializer.force_instantiate();
   }
};

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

} // namespace detail

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   digamma(T x, const Policy&)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   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_<114> >
      >,
      mpl::int_<0>,
      typename mpl::if_<
         mpl::less<precision_type, mpl::int_<25> >,
         mpl::int_<24>,
         typename mpl::if_<
            mpl::less<precision_type, mpl::int_<54> >,
            mpl::int_<53>,
            typename mpl::if_<
               mpl::less<precision_type, mpl::int_<65> >,
               mpl::int_<64>,
               mpl::int_<113>
            >::type
         >::type
      >::type
   >::type tag_type;

   typedef typename policies::normalise<
      Policy,
      policies::promote_float<false>,
      policies::promote_double<false>,
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   // Force initialization of constants:
   detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate();

   return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
      static_cast<value_type>(x),
      static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::digamma<%1%>(%1%)");
}

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

} // namespace math
} // namespace boost

#ifdef _MSC_VER
#pragma warning(pop)
#endif

#endif


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  Copyright (c) 2006 John Maddock
//  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)
//
//  History:
//  XZ wrote the original of this file as part of the Google
//  Summer of Code 2006.  JM modified it to fit into the
//  Boost.Math conceptual framework better, and to ensure
//  that the code continues to work no matter how many digits
//  type T has.

#ifndef BOOST_MATH_ELLINT_1_HPP
#define BOOST_MATH_ELLINT_1_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/ellint_rf.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/workaround.hpp>
#include <boost/math/special_functions/round.hpp>

// Elliptic integrals (complete and incomplete) of the first kind
// Carlson, Numerische Mathematik, vol 33, 1 (1979)

namespace boost { namespace math {

template <class T1, class T2, class Policy>
typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol);

namespace detail{

template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol);

// Elliptic integral (Legendre form) of the first kind
template <typename T, typename Policy>
T ellint_f_imp(T phi, T k, const Policy& pol)
{
    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    static const char* function = "boost::math::ellint_f<%1%>(%1%,%1%)";
    BOOST_MATH_INSTRUMENT_VARIABLE(phi);
    BOOST_MATH_INSTRUMENT_VARIABLE(k);
    BOOST_MATH_INSTRUMENT_VARIABLE(function);

    if (abs(k) > 1)
    {
       return policies::raise_domain_error<T>(function,
            "Got k = %1%, function requires |k| <= 1", k, pol);
    }

    bool invert = false;
    if(phi < 0)
    {
       BOOST_MATH_INSTRUMENT_VARIABLE(phi);
       phi = fabs(phi);
       invert = true;
    }

    T result;

    if(phi >= tools::max_value<T>())
    {
       // Need to handle infinity as a special case:
       result = policies::raise_overflow_error<T>(function, 0, pol);
       BOOST_MATH_INSTRUMENT_VARIABLE(result);
    }
    else if(phi > 1 / tools::epsilon<T>())
    {
       // Phi is so large that phi%pi is necessarily zero (or garbage),
       // just return the second part of the duplication formula:
       result = 2 * phi * ellint_k_imp(k, pol) / constants::pi<T>();
       BOOST_MATH_INSTRUMENT_VARIABLE(result);
    }
    else
    {
       // Carlson's algorithm works only for |phi| <= pi/2,
       // use the integrand's periodicity to normalize phi
       //
       // Xiaogang's original code used a cast to long long here
       // but that fails if T has more digits than a long long,
       // so rewritten to use fmod instead:
       //
       BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
       T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
       BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
       T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
       BOOST_MATH_INSTRUMENT_VARIABLE(m);
       int s = 1;
       if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
       {
          m += 1;
          s = -1;
          rphi = constants::half_pi<T>() - rphi;
          BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
       }
       T sinp = sin(rphi);
       sinp *= sinp;
       T cosp = cos(rphi);
       cosp *= cosp;
       BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
       BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
       if(sinp > tools::min_value<T>())
       {
          //
          // Use http://dlmf.nist.gov/19.25#E5, note that
          // c-1 simplifies to cot^2(rphi) which avoid cancellation:
          //
          T c = 1 / sinp;
          result = rphi == 0 ? static_cast<T>(0) : static_cast<T>(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol));
       }
       else
          result = s * sin(rphi);
       BOOST_MATH_INSTRUMENT_VARIABLE(result);
       if(m != 0)
       {
          result += m * ellint_k_imp(k, pol);
          BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
    }
    return invert ? T(-result) : result;
}

// Complete elliptic integral (Legendre form) of the first kind
template <typename T, typename Policy>
T ellint_k_imp(T k, const Policy& pol)
{
    BOOST_MATH_STD_USING
    using namespace boost::math::tools;

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

    if (abs(k) > 1)
    {
       return policies::raise_domain_error<T>(function,
            "Got k = %1%, function requires |k| <= 1", k, pol);
    }
    if (abs(k) == 1)
    {
       return policies::raise_overflow_error<T>(function, 0, pol);
    }

    T x = 0;
    T y = 1 - k * k;
    T z = 1;
    T value = ellint_rf_imp(x, y, z, pol);

    return value;
}

template <typename T, typename Policy>
inline typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, const mpl::true_&)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_k_imp(static_cast<value_type>(k), pol), "boost::math::ellint_1<%1%>(%1%)");
}

template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const mpl::false_&)
{
   return boost::math::ellint_1(k, phi, policies::policy<>());
}

}

// Complete elliptic integral (Legendre form) of the first kind
template <typename T>
inline typename tools::promote_args<T>::type ellint_1(T k)
{
   return ellint_1(k, policies::policy<>());
}

// Elliptic integral (Legendre form) of the first kind
template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol)
{
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_f_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_1<%1%>(%1%,%1%)");
}

template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi)
{
   typedef typename policies::is_policy<T2>::type tag_type;
   return detail::ellint_1(k, phi, tag_type());
}

}} // namespaces

#endif // BOOST_MATH_ELLINT_1_HPP


code:
//  Copyright (c) 2006 Xiaogang Zhang
//  Copyright (c) 2006 John Maddock
//  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)
//
//  History:
//  XZ wrote the original of this file as part of the Google
//  Summer of Code 2006.  JM modified it to fit into the
//  Boost.Math conceptual framework better, and to ensure
//  that the code continues to work no matter how many digits
//  type T has.

#ifndef BOOST_MATH_ELLINT_2_HPP
#define BOOST_MATH_ELLINT_2_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/ellint_rf.hpp>
#include <boost/math/special_functions/ellint_rd.hpp>
#include <boost/math/special_functions/ellint_rg.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/tools/workaround.hpp>
#include <boost/math/special_functions/round.hpp>

// Elliptic integrals (complete and incomplete) of the second kind
// Carlson, Numerische Mathematik, vol 33, 1 (1979)

namespace boost { namespace math {

template <class T1, class T2, class Policy>
typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol);

namespace detail{

template <typename T, typename Policy>
T ellint_e_imp(T k, const Policy& pol);

// Elliptic integral (Legendre form) of the second kind
template <typename T, typename Policy>
T ellint_e_imp(T phi, T k, const Policy& pol)
{
    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    bool invert = false;
    if(phi < 0)
    {
       phi = fabs(phi);
       invert = true;
    }

    T result;

    if(phi >= tools::max_value<T>())
    {
       // Need to handle infinity as a special case:
       result = policies::raise_overflow_error<T>("boost::math::ellint_e<%1%>(%1%,%1%)", 0, pol);
    }
    else if(phi > 1 / tools::epsilon<T>())
    {
       // Phi is so large that phi%pi is necessarily zero (or garbage),
       // just return the second part of the duplication formula:
       result = 2 * phi * ellint_e_imp(k, pol) / constants::pi<T>();
    }
    else if(k == 0)
    {
       return invert ? T(-phi) : phi;
    }
    else if(fabs(k) == 1)
    {
       return invert ? T(-sin(phi)) : T(sin(phi));
    }
    else
    {
       // Carlson's algorithm works only for |phi| <= pi/2,
       // use the integrand's periodicity to normalize phi
       //
       // Xiaogang's original code used a cast to long long here
       // but that fails if T has more digits than a long long,
       // so rewritten to use fmod instead:
       //
       T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
       T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
       int s = 1;
       if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
       {
          m += 1;
          s = -1;
          rphi = constants::half_pi<T>() - rphi;
       }
       T k2 = k * k;
       if(k2 > 1)
       {
          return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
       }
       else if(rphi < tools::root_epsilon<T>())
       {
          // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/
          result = s * rphi;
       }
       else
       {
          // http://dlmf.nist.gov/19.25#E10
          T sinp = sin(rphi);
          T cosp = cos(rphi);
          T c = 1 / (sinp * sinp);
          T cm1 = cosp * cosp / (sinp * sinp);  // c - 1
          result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2))));
       }
       if(m != 0)
          result += m * ellint_e_imp(k, pol);
    }
    return invert ? T(-result) : result;
}

// Complete elliptic integral (Legendre form) of the second kind
template <typename T, typename Policy>
T ellint_e_imp(T k, const Policy& pol)
{
    BOOST_MATH_STD_USING
    using namespace boost::math::tools;

    if (abs(k) > 1)
    {
       return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)",
            "Got k = %1%, function requires |k| <= 1", k, pol);
    }
    if (abs(k) == 1)
    {
        return static_cast<T>(1);
    }

    T x = 0;
    T t = k * k;
    T y = 1 - t;
    T z = 1;
    T value = 2 * ellint_rg_imp(x, y, z, pol);

    return value;
}

template <typename T, typename Policy>
inline typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const mpl::true_&)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_e_imp(static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%)");
}

// Elliptic integral (Legendre form) of the second kind
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const mpl::false_&)
{
   return boost::math::ellint_2(k, phi, policies::policy<>());
}

} // detail

// Complete elliptic integral (Legendre form) of the second kind
template <typename T>
inline typename tools::promote_args<T>::type ellint_2(T k)
{
   return ellint_2(k, policies::policy<>());
}

// Elliptic integral (Legendre form) of the second kind
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)
{
   typedef typename policies::is_policy<T2>::type tag_type;
   return detail::ellint_2(k, phi, tag_type());
}

template <class T1, class T2, class Policy>
inline typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol)
{
   typedef typename tools::promote_args<T1, T2>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_e_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%,%1%)");
}

}} // namespaces

#endif // BOOST_MATH_ELLINT_2_HPP