X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=boost%2Fmath%2Fspecial_functions%2Fgamma.hpp;fp=boost%2Fmath%2Fspecial_functions%2Fgamma.hpp;h=8e661a1fdfbfd1f3d484326ab9ef1a9767eedede;hp=f809fa6f40634bd346b8275bfd815e6e9f45661c;hb=2d71eb92104693ca9baa5a2e1c23eeca776d8fd3;hpb=da57529b92adbb7ae74a89861cb39fb35ac7c62d

diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index f809fa6..8e661a1 100644
--- a/boost/math/special_functions/gamma.hpp
+++ b/boost/math/special_functions/gamma.hpp
@@ -1,3 +1,4 @@
+
 //  Copyright John Maddock 2006-7.
 //  Copyright Paul A. Bristow 2007.
 
@@ -13,17 +14,6 @@
 #endif
 
 #include <boost/config.hpp>
-#ifdef BOOST_MSVC
-# pragma warning(push)
-# pragma warning(disable: 4127 4701)
-//  // For lexical_cast, until fixed in 1.35?
-//  // conditional expression is constant &
-//  // Potentially uninitialized local variable 'name' used
-#endif
-#include <boost/lexical_cast.hpp>
-#ifdef BOOST_MSVC
-# pragma warning(pop)
-#endif
 #include <boost/math/tools/series.hpp>
 #include <boost/math/tools/fraction.hpp>
 #include <boost/math/tools/precision.hpp>
@@ -49,12 +39,6 @@
 #include <boost/config/no_tr1/cmath.hpp>
 #include <algorithm>
 
-#ifdef BOOST_MATH_INSTRUMENT
-#include <iostream>
-#include <iomanip>
-#include <typeinfo>
-#endif
-
 #ifdef BOOST_MSVC
 # pragma warning(push)
 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
@@ -126,8 +110,8 @@ T sinpx(T z)
 //
 // tgamma(z), with Lanczos support:
 //
-template <class T, class Policy, class L>
-T gamma_imp(T z, const Policy& pol, const L& l)
+template <class T, class Policy, class Lanczos>
+T gamma_imp(T z, const Policy& pol, const Lanczos& l)
 {
    BOOST_MATH_STD_USING
 
@@ -150,6 +134,7 @@ T gamma_imp(T z, const Policy& pol, const L& l)
       if(z <= -20)
       {
          result = gamma_imp(T(-z), pol, l) * sinpx(z);
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
          if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
             return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
          result = -boost::math::constants::pi<T>() / result;
@@ -157,6 +142,7 @@ T gamma_imp(T z, const Policy& pol, const L& l)
             return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
          if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
             return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
          return result;
       }
 
@@ -167,29 +153,41 @@ T gamma_imp(T z, const Policy& pol, const L& l)
          z += 1;
       }
    }
+   BOOST_MATH_INSTRUMENT_VARIABLE(result);
    if((floor(z) == z) && (z < max_factorial<T>::value))
    {
       result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
+      BOOST_MATH_INSTRUMENT_VARIABLE(result);
    }
    else
    {
-      result *= L::lanczos_sum(z);
-      if(z * log(z) > tools::log_max_value<T>())
+      result *= Lanczos::lanczos_sum(z);
+      T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
+      T lzgh = log(zgh);
+      BOOST_MATH_INSTRUMENT_VARIABLE(result);
+      BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
+      if(z * lzgh > tools::log_max_value<T>())
       {
          // we're going to overflow unless this is done with care:
-         T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
-         if(log(zgh) * z / 2 > tools::log_max_value<T>())
+         BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
+         if(lzgh * z / 2 > tools::log_max_value<T>())
             return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
          T hp = pow(zgh, (z / 2) - T(0.25));
+         BOOST_MATH_INSTRUMENT_VARIABLE(hp);
          result *= hp / exp(zgh);
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
          if(tools::max_value<T>() / hp < result)
             return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
          result *= hp;
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
       else
       {
-         T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
+         BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
+         BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
+         BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
          result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
    }
    return result;
@@ -197,8 +195,8 @@ T gamma_imp(T z, const Policy& pol, const L& l)
 //
 // lgamma(z) with Lanczos support:
 //
-template <class T, class Policy, class L>
-T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
+template <class T, class Policy, class Lanczos>
+T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
 {
 #ifdef BOOST_MATH_INSTRUMENT
    static bool b = false;
@@ -249,9 +247,9 @@ T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
             >,
             mpl::int_<113>, mpl::int_<0> >::type
           >::type tag_type;
-      result = lgamma_small_imp<T>(z, z - 1, z - 2, tag_type(), pol, l);
+      result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
    }
-   else if((z >= 3) && (z < 100))
+   else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
    {
       // taking the log of tgamma reduces the error, no danger of overflow here:
       result = log(gamma_imp(z, pol, l));
@@ -259,10 +257,10 @@ T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
    else
    {
       // regular evaluation:
-      T zgh = static_cast<T>(z + L::g() - boost::math::constants::half<T>());
+      T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
       result = log(zgh) - 1;
       result *= z - 0.5f;
-      result += log(L::lanczos_sum_expG_scaled(z));
+      result += log(Lanczos::lanczos_sum_expG_scaled(z));
    }
 
    if(sign)
@@ -333,7 +331,7 @@ inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
    T factor = policies::get_epsilon<T, Policy>();
    T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
-   policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
+   policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
    return result;
 }
 
@@ -350,7 +348,7 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
       return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
    if(z <= -20)
    {
-      T result = gamma_imp(-z, pol, l) * sinpx(z);
+      T result = gamma_imp(T(-z), pol, l) * sinpx(z);
       if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
          return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
       result = -boost::math::constants::pi<T>() / result;
@@ -418,7 +416,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*si
    }
    else if((z != 1) && (z != 2))
    {
-      T limit = (std::max)(z+1, T(10));
+      T limit = (std::max)(T(z+1), T(10));
       T prefix = z * log(limit) - limit;
       T sum = detail::lower_gamma_series(z, limit, pol) / z;
       sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
@@ -432,8 +430,8 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*si
 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
 // used by the upper incomplete gamma with z < 1:
 //
-template <class T, class Policy, class L>
-T tgammap1m1_imp(T dz, Policy const& pol, const L& l)
+template <class T, class Policy, class Lanczos>
+T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
 {
    BOOST_MATH_STD_USING
 
@@ -445,7 +443,7 @@ T tgammap1m1_imp(T dz, Policy const& pol, const L& l)
          mpl::greater<precision_type, mpl::int_<113> >
       >,
       typename mpl::if_<
-         is_same<L, lanczos::lanczos24m113>,
+         is_same<Lanczos, lanczos::lanczos24m113>,
          mpl::int_<113>,
          mpl::int_<0>
       >::type,
@@ -500,7 +498,7 @@ inline T tgammap1m1_imp(T dz, Policy const& pol,
    // algebra isn't easy for the general case....
    // Start by subracting 1 from tgamma:
    //
-   T result = gamma_imp(1 + dz, pol, l) - 1;
+   T result = gamma_imp(T(1 + dz), pol, l) - 1;
    BOOST_MATH_INSTRUMENT_CODE(result);
    //
    // Test the level of cancellation error observed: we loose one bit
@@ -597,13 +595,13 @@ T full_igamma_prefix(T a, T z, const Policy& pol)
 // Compute (z^a)(e^-z)/tgamma(a)
 // most if the error occurs in this function:
 //
-template <class T, class Policy, class L>
-T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
+template <class T, class Policy, class Lanczos>
+T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
 {
    BOOST_MATH_STD_USING
-   T agh = a + static_cast<T>(L::g()) - T(0.5);
+   T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
    T prefix;
-   T d = ((z - a) - static_cast<T>(L::g()) + T(0.5)) / agh;
+   T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
 
    if(a < 1)
    {
@@ -631,7 +629,7 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
    else if((fabs(d*d*a) <= 100) && (a > 150))
    {
       // special case for large a and a ~ z.
-      prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - L::g()) / agh;
+      prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
       prefix = exp(prefix);
    }
    else
@@ -673,7 +671,7 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
          prefix = pow(z / agh, a) * exp(amz);
       }
    }
-   prefix *= sqrt(agh / boost::math::constants::e<T>()) / L::lanczos_sum_expG_scaled(a);
+   prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
    return prefix;
 }
 //
@@ -748,7 +746,7 @@ inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool
       *pderivative = p / (*pgam * exp(x));
    T init_value = invert ? *pgam : 0;
    result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
-   policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
+   policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
    if(invert)
       result = -result;
    return result;
@@ -835,12 +833,12 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
 
    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
 
-   T result;
+   T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
 
    BOOST_ASSERT((p_derivative == 0) || (normalised == true));
 
    bool is_int, is_half_int;
-   bool is_small_a = (a < 30) && (a <= x + 1);
+   bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
    if(is_small_a)
    {
       T fa = floor(a);
@@ -1077,11 +1075,11 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
 //
 // Ratios of two gamma functions:
 //
-template <class T, class Policy, class L>
-T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const L&)
+template <class T, class Policy, class Lanczos>
+T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&)
 {
    BOOST_MATH_STD_USING
-   T zgh = z + L::g() - constants::half<T>();
+   T zgh = z + Lanczos::g() - constants::half<T>();
    T result;
    if(fabs(delta) < 10)
    {
@@ -1092,7 +1090,7 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const L&)
       result = pow(zgh / (zgh + delta), z - constants::half<T>());
    }
    result *= pow(constants::e<T>() / (zgh + delta), delta);
-   result *= L::lanczos_sum(z) / L::lanczos_sum(z + delta);
+   result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
    return result;
 }
 //
@@ -1192,6 +1190,68 @@ T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
    return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
 }
 
+template <class T, class Policy>
+T tgamma_ratio_imp(T x, T y, const Policy& pol)
+{
+   BOOST_MATH_STD_USING
+
+   if((x <= tools::min_value<T>()) || (boost::math::isinf)(x))
+      policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
+   if((y <= tools::min_value<T>()) || (boost::math::isinf)(y))
+      policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);
+
+   if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
+   {
+      // Rather than subtracting values, lets just call the gamma functions directly:
+      return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+   }
+   T prefix = 1;
+   if(x < 1)
+   {
+      if(y < 2 * max_factorial<T>::value)
+      {
+         // We need to sidestep on x as well, otherwise we'll underflow
+         // before we get to factor in the prefix term:
+         prefix /= x;
+         x += 1;
+         while(y >=  max_factorial<T>::value)
+         {
+            y -= 1;
+            prefix /= y;
+         }
+         return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+      }
+      //
+      // result is almost certainly going to underflow to zero, try logs just in case:
+      //
+      return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+   }
+   if(y < 1)
+   {
+      if(x < 2 * max_factorial<T>::value)
+      {
+         // We need to sidestep on y as well, otherwise we'll overflow
+         // before we get to factor in the prefix term:
+         prefix *= y;
+         y += 1;
+         while(x >= max_factorial<T>::value)
+         {
+            x -= 1;
+            prefix *= x;
+         }
+         return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
+      }
+      //
+      // Result will almost certainly overflow, try logs just in case:
+      //
+      return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
+   }
+   //
+   // Regular case, x and y both large and similar in magnitude:
+   //
+   return boost::math::tgamma_delta_ratio(x, y - x, pol);
+}
+
 template <class T, class Policy>
 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
 {
@@ -1243,6 +1303,101 @@ inline typename tools::promote_args<T>::type
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)");
 }
 
+template <class T, class Policy>
+struct igamma_initializer
+{
+   struct init
+   {
+      init()
+      {
+         typedef typename policies::precision<T, Policy>::type precision_type;
+
+         typedef typename mpl::if_<
+            mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
+            mpl::greater<precision_type, mpl::int_<113> > >,
+            mpl::int_<0>,
+            typename mpl::if_<
+               mpl::less_equal<precision_type, mpl::int_<53> >,
+               mpl::int_<53>,
+               typename mpl::if_<
+                  mpl::less_equal<precision_type, mpl::int_<64> >,
+                  mpl::int_<64>,
+                  mpl::int_<113>
+               >::type
+            >::type
+         >::type tag_type;
+
+         do_init(tag_type());
+      }
+      template <int N>
+      static void do_init(const mpl::int_<N>&)
+      {
+         boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
+      }
+      static void do_init(const mpl::int_<53>&){}
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
+
+template <class T, class Policy>
+struct lgamma_initializer
+{
+   struct init
+   {
+      init()
+      {
+         typedef typename policies::precision<T, Policy>::type precision_type;
+         typedef typename mpl::if_<
+            mpl::and_<
+               mpl::less_equal<precision_type, mpl::int_<64> >, 
+               mpl::greater<precision_type, mpl::int_<0> > 
+            >,
+            mpl::int_<64>,
+            typename mpl::if_<
+               mpl::and_<
+                  mpl::less_equal<precision_type, mpl::int_<113> >,
+                  mpl::greater<precision_type, mpl::int_<0> > 
+               >,
+               mpl::int_<113>, mpl::int_<0> >::type
+             >::type tag_type;
+         do_init(tag_type());
+      }
+      static void do_init(const mpl::int_<64>&)
+      {
+         boost::math::lgamma(static_cast<T>(2.5), Policy());
+         boost::math::lgamma(static_cast<T>(1.25), Policy());
+         boost::math::lgamma(static_cast<T>(1.75), Policy());
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         boost::math::lgamma(static_cast<T>(2.5), Policy());
+         boost::math::lgamma(static_cast<T>(1.25), Policy());
+         boost::math::lgamma(static_cast<T>(1.5), Policy());
+         boost::math::lgamma(static_cast<T>(1.75), Policy());
+      }
+      static void do_init(const mpl::int_<0>&)
+      {
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy>
+const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
+
 template <class T1, class T2, class Policy>
 inline typename tools::promote_args<T1, T2>::type
    tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
@@ -1250,13 +1405,16 @@ inline typename tools::promote_args<T1, T2>::type
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename tools::promote_args<T1, T2>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
-   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
    typedef typename policies::normalise<
       Policy, 
       policies::promote_float<false>, 
       policies::promote_double<false>, 
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
+
+   igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), false, true,
@@ -1270,6 +1428,7 @@ inline typename tools::promote_args<T1, T2>::type
    return tgamma(a, z, policies::policy<>(), tag);
 }
 
+
 } // namespace detail
 
 template <class T>
@@ -1293,6 +1452,9 @@ inline typename tools::promote_args<T>::type
       policies::promote_double<false>, 
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
+
+   detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)");
 }
 
@@ -1372,7 +1534,7 @@ inline typename tools::promote_args<T1, T2>::type
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename tools::promote_args<T1, T2>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
-   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
    typedef typename policies::normalise<
       Policy, 
       policies::promote_float<false>, 
@@ -1380,6 +1542,8 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), false, false,
@@ -1401,7 +1565,7 @@ inline typename tools::promote_args<T1, T2>::type
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename tools::promote_args<T1, T2>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
-   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
    typedef typename policies::normalise<
       Policy, 
       policies::promote_float<false>, 
@@ -1409,6 +1573,8 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), true, true,
@@ -1430,7 +1596,7 @@ inline typename tools::promote_args<T1, T2>::type
    BOOST_FPU_EXCEPTION_GUARD
    typedef typename tools::promote_args<T1, T2>::type result_type;
    typedef typename policies::evaluation<result_type, Policy>::type value_type;
-   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+   // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
    typedef typename policies::normalise<
       Policy, 
       policies::promote_float<false>, 
@@ -1438,6 +1604,8 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
+   detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
       detail::gamma_incomplete_imp(static_cast<value_type>(a),
       static_cast<value_type>(z), true, false,
@@ -1486,7 +1654,7 @@ inline typename tools::promote_args<T1, T2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
 }
 template <class T1, class T2>
 inline typename tools::promote_args<T1, T2>::type