view release on metacpan or  search on metacpan

src/boost/math/special_functions/gamma.hpp  view on Meta::CPAN

   {
      result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
      BOOST_MATH_INSTRUMENT_VARIABLE(result);
   }
   else
   {
      result *= Lanczos::lanczos_sum(z);
      BOOST_MATH_INSTRUMENT_VARIABLE(result);
      BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
      if(z * log(z) > tools::log_max_value<T>())
      {
         // we're going to overflow unless this is done with care:
         T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
         BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
         if(log(zgh) * 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>(Lanczos::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;
}
//
// lgamma(z) with Lanczos support:
//
template <class T, class Policy, class Lanczos>
T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
{
#ifdef BOOST_MATH_INSTRUMENT
   static bool b = false;
   if(!b)
   {
      std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
      b = true;
   }
#endif

   BOOST_MATH_STD_USING

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

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

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

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

//
// Incomplete gamma functions follow:
//
template <class T>
struct upper_incomplete_gamma_fract
{
private:
   T z, a;