view release on metacpan or  search on metacpan

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

         BOOST_MATH_BIG_CONSTANT(T, 113, -1077042281708.42654526404581272546244),
         BOOST_MATH_BIG_CONSTANT(T, 113, -68028222642.1941480871395695677675137)
      };
      static const T Q[20] = {    
         BOOST_MATH_BIG_CONSTANT(T, 113, 1),
         BOOST_MATH_BIG_CONSTANT(T, 113, 168.542326331163836642960118190147311),
         BOOST_MATH_BIG_CONSTANT(T, 113, 12535.7237814586576783518249115343619),
         BOOST_MATH_BIG_CONSTANT(T, 113, 544891.263372016404143120911148640627),
         BOOST_MATH_BIG_CONSTANT(T, 113, 15454474.7241010258634446523045237762),
         BOOST_MATH_BIG_CONSTANT(T, 113, 302495899.896629522673410325891717381),
         BOOST_MATH_BIG_CONSTANT(T, 113, 4215565948.38886507646911672693270307),
         BOOST_MATH_BIG_CONSTANT(T, 113, 42552409471.7951815668506556705733344),
         BOOST_MATH_BIG_CONSTANT(T, 113, 313592377066.753173979584098301610186),
         BOOST_MATH_BIG_CONSTANT(T, 113, 1688763640223.4541980740597514904542),
         BOOST_MATH_BIG_CONSTANT(T, 113, 6610992294901.59589748057620192145704),
         BOOST_MATH_BIG_CONSTANT(T, 113, 18601637235659.6059890851321772682606),
         BOOST_MATH_BIG_CONSTANT(T, 113, 36944278231087.2571020964163402941583),
         BOOST_MATH_BIG_CONSTANT(T, 113, 50425858518481.7497071917028793820058),
         BOOST_MATH_BIG_CONSTANT(T, 113, 45508060902865.0899967797848815980644),
         BOOST_MATH_BIG_CONSTANT(T, 113, 25649955002765.3817331501988304758142),
         BOOST_MATH_BIG_CONSTANT(T, 113, 8259575619094.6518520988612711292331),
         BOOST_MATH_BIG_CONSTANT(T, 113, 1299981487496.12607474362723586264515),
         BOOST_MATH_BIG_CONSTANT(T, 113, 70242279152.8241187845178443118302693),
         BOOST_MATH_BIG_CONSTANT(T, 113, -37633302.9409263839042721539363416685)
      };
      T recip = 1 / z;
      result = 1 + tools::evaluate_polynomial(P, recip)
         / tools::evaluate_polynomial(Q, recip);
      result *= exp(-z) * recip;
   }
   else
   {
      result = 0;
   }
   return result;
}

template <class T>
struct expint_fraction
{
   typedef std::pair<T,T> result_type;
   expint_fraction(unsigned n_, T z_) : b(n_ + z_), i(-1), n(n_){}
   std::pair<T,T> operator()()
   {
      std::pair<T,T> result = std::make_pair(-static_cast<T>((i+1) * (n+i)), b);
      b += 2;
      ++i;
      return result;
   }
private:
   T b;
   int i;
   unsigned n;
};

template <class T, class Policy>
inline T expint_as_fraction(unsigned n, T z, const Policy& pol)
{
   BOOST_MATH_STD_USING
   BOOST_MATH_INSTRUMENT_VARIABLE(z)
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
   expint_fraction<T> f(n, z);
   T result = tools::continued_fraction_b(
      f, 
      boost::math::policies::get_epsilon<T, Policy>(),
      max_iter);
   policies::check_series_iterations<T>("boost::math::expint_continued_fraction<%1%>(unsigned,%1%)", max_iter, pol);
   BOOST_MATH_INSTRUMENT_VARIABLE(result)
   BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
   result = exp(-z) / result;
   BOOST_MATH_INSTRUMENT_VARIABLE(result)
   return result;
}

template <class T>
struct expint_series
{
   typedef T result_type;
   expint_series(unsigned k_, T z_, T x_k_, T denom_, T fact_) 
      : k(k_), z(z_), x_k(x_k_), denom(denom_), fact(fact_){}
   T operator()()
   {
      x_k *= -z;
      denom += 1;
      fact *= ++k;
      return x_k / (denom * fact);
   }
private:
   unsigned k;
   T z;
   T x_k;
   T denom;
   T fact;
};

template <class T, class Policy>
inline T expint_as_series(unsigned n, T z, const Policy& pol)
{
   BOOST_MATH_STD_USING
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();

   BOOST_MATH_INSTRUMENT_VARIABLE(z)

   T result = 0;
   T x_k = -1;
   T denom = T(1) - n;
   T fact = 1;
   unsigned k = 0;
   for(; k < n - 1;)
   {
      result += x_k / (denom * fact);
      denom += 1;
      x_k *= -z;
      fact *= ++k;
   }
   BOOST_MATH_INSTRUMENT_VARIABLE(result)
   result += pow(-z, static_cast<T>(n - 1)) 
      * (boost::math::digamma(static_cast<T>(n)) - log(z)) / fact;
   BOOST_MATH_INSTRUMENT_VARIABLE(result)

   expint_series<T> s(k, z, x_k, denom, fact);
   result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
   policies::check_series_iterations<T>("boost::math::expint_series<%1%>(unsigned,%1%)", max_iter, pol);
   BOOST_MATH_INSTRUMENT_VARIABLE(result)
   BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
   return result;
}

template <class T, class Policy, class Tag>
T expint_imp(unsigned n, T z, const Policy& pol, const Tag& tag)
{
   BOOST_MATH_STD_USING
   static const char* function = "boost::math::expint<%1%>(unsigned, %1%)";
   if(z < 0)
      return policies::raise_domain_error<T>(function, "Function requires z >= 0 but got %1%.", z, pol);
   if(z == 0)
      return n == 1 ? policies::raise_overflow_error<T>(function, 0, pol) : T(1 / (static_cast<T>(n - 1)));

   T result;

   bool f;
   if(n < 3)
   {
      f = z < 0.5;
   }
   else
   {
      f = z < (static_cast<T>(n - 2) / static_cast<T>(n - 1));
   }
#ifdef BOOST_MSVC
#  pragma warning(push)
#  pragma warning(disable:4127) // conditional expression is constant
#endif
   if(n == 0)
      result = exp(-z) / z;
   else if((n == 1) && (Tag::value))
   {
      result = expint_1_rational(z, tag);
   }
   else if(f)
      result = expint_as_series(n, z, pol);
   else
      result = expint_as_fraction(n, z, pol);
#ifdef BOOST_MSVC
#  pragma warning(pop)
#endif

   return result;
}

template <class T>
struct expint_i_series
{
   typedef T result_type;
   expint_i_series(T z_) : k(0), z_k(1), z(z_){}
   T operator()()
   {
      z_k *= z / ++k;
      return z_k / k;
   }
private:
   unsigned k;
   T z_k;
   T z;
};

template <class T, class Policy>
T expint_i_as_series(T z, const Policy& pol)
{
   BOOST_MATH_STD_USING
   T result = log(z); // (log(z) - log(1 / z)) / 2;
   result += constants::euler<T>();
   expint_i_series<T> s(z);
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
   result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
   policies::check_series_iterations<T>("boost::math::expint_i_series<%1%>(%1%)", max_iter, pol);
   return result;
}

template <class T, class Policy, class Tag>
T expint_i_imp(T z, const Policy& pol, const Tag& tag)
{
   static const char* function = "boost::math::expint<%1%>(%1%)";
   if(z < 0)
      return -expint_imp(1, T(-z), pol, tag);
   if(z == 0)
      return -policies::raise_overflow_error<T>(function, 0, pol);
   return expint_i_as_series(z, pol);
}

template <class T, class Policy>
T expint_i_imp(T z, const Policy& pol, const mpl::int_<53>& tag)
{
   BOOST_MATH_STD_USING
   static const char* function = "boost::math::expint<%1%>(%1%)";
   if(z < 0)
      return -expint_imp(1, T(-z), pol, tag);
   if(z == 0)
      return -policies::raise_overflow_error<T>(function, 0, pol);

   T result;

   if(z <= 6)
   {
      // Maximum Deviation Found:                     2.852e-18
      // Expected Error Term:                         2.852e-18
      // Max Error found at double precision =        Poly: 2.636335e-16   Cheb: 4.187027e-16
      static const T P[10] = {    
         BOOST_MATH_BIG_CONSTANT(T, 53, 2.98677224343598593013),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.356343618769377415068),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.780836076283730801839),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.114670926327032002811),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.0499434773576515260534),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.00726224593341228159561),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.00115478237227804306827),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.000116419523609765200999),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.798296365679269702435e-5),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.2777056254402008721e-6)
      };
      static const T Q[8] = {    
         BOOST_MATH_BIG_CONSTANT(T, 53, 1),
         BOOST_MATH_BIG_CONSTANT(T, 53, -1.17090412365413911947),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.62215109846016746276),
         BOOST_MATH_BIG_CONSTANT(T, 53, -0.195114782069495403315),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.0391523431392967238166),
         BOOST_MATH_BIG_CONSTANT(T, 53, -0.00504800158663705747345),
         BOOST_MATH_BIG_CONSTANT(T, 53, 0.000389034007436065401822),
         BOOST_MATH_BIG_CONSTANT(T, 53, -0.138972589601781706598e-4)
      };

      static const T c1 = BOOST_MATH_BIG_CONSTANT(T, 53, 1677624236387711.0);
      static const T c2 = BOOST_MATH_BIG_CONSTANT(T, 53, 4503599627370496.0);
      static const T r1 = static_cast<T>(c1 / c2);
      static const T r2 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.131401834143860282009280387409357165515556574352422001206362e-16);
      static const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 53, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));
      T t = (z / 3) - 1;