expint.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 1,515 行 · 第 1/4 页

{
   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::digits<T, Policy>(), max_iter, result);
   policies::check_series_iterations("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) : 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::digits<T, Policy>(), max_iter, result);
   policies::check_series_iterations("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, -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, -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] = {    
         2.98677224343598593013L,
         0.356343618769377415068L,
         0.780836076283730801839L,
         0.114670926327032002811L,
         0.0499434773576515260534L,
         0.00726224593341228159561L,
         0.00115478237227804306827L,
         0.000116419523609765200999L,
         0.798296365679269702435e-5L,
         0.2777056254402008721e-6L
      };
      static const T Q[8] = {    
         1L,
         -1.17090412365413911947L,
         0.62215109846016746276L,
         -0.195114782069495403315L,
         0.0391523431392967238166L,
         -0.00504800158663705747345L,
         0.000389034007436065401822L,
         -0.138972589601781706598e-4L
      };

      static const T r1 = static_cast<T>(1677624236387711.0L) / 4503599627370496uLL;
      static const T r2 = 0.131401834143860282009280387409357165515556574352422001206362e-16L;
      static const T r = static_cast<T>(0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392L);
      T t = (z / 3) - 1;
      result = tools::evaluate_polynomial(P, t) 
         / tools::evaluate_polynomial(Q, t);
      t = (z - r1) - r2;
      result *= t;
      if(fabs(t) < 0.1)
      {
         result += boost::math::log1p(t / r);
      }
      else
      {
         result += log(z / r);
      }
   }
   else if (z <= 10)
   {
      // Maximum Deviation Found:                     6.546e-17
      // Expected Error Term:                         6.546e-17
      // Max Error found at double precision =        Poly: 6.890169e-17   Cheb: 6.772128e-17
      static const T Y = 1.158985137939453125F;
      static const T P[8] = {    
         0.00139324086199402804173L,
         -0.0349921221823888744966L,
         -0.0264095520754134848538L,
         -0.00761224003005476438412L,
         -0.00247496209592143627977L,
         -0.000374885917942100256775L,
         -0.554086272024881826253e-4L,
         -0.396487648924804510056e-5L
      };
      static const T Q[8] = {    
         1L,
         0.744625566823272107711L,
         0.329061095011767059236L,
         0.100128624977313872323L,
         0.0223851099128506347278L,
         0.00365334190742316650106L,
         0.000402453408512476836472L,
         0.263649630720255691787e-4L
      };
      T t = z / 2 - 4;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      result *= exp(z) / z;
      result += z;
   }
   else if(z <= 20)
   {
      // Maximum Deviation Found:                     1.843e-17
      // Expected Error Term:                         -1.842e-17
      // Max Error found at double precision =        Poly: 4.375868e-17   Cheb: 5.860967e-17

      static const T Y = 1.0869731903076171875F;
      static const T P[9] = {    
         -0.00893891094356945667451L,
         -0.0484607730127134045806L,
         -0.0652810444222236895772L,
         -0.0478447572647309671455L,
         -0.0226059218923777094596L,
         -0.00720603636917482065907L,
         -0.00155941947035972031334L,
         -0.000209750022660200888349L,
         -0.138652200349182596186e-4L
      };
      static const T Q[9] = {    
         1L,
         1.97017214039061194971L,
         1.86232465043073157508L,
         1.09601437090337519977L,
         0.438873285773088870812L,
         0.122537731979686102756L,
         0.0233458478275769288159L,
         0.00278170769163303669021L,
         0.000159150281166108755531L
      };
      T t = z / 5 - 3;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      result *= exp(z) / z;
      result += z;
   }
   else if(z <= 40)
   {
      // Maximum Deviation Found:                     5.102e-18
      // Expected Error Term:                         5.101e-18
      // Max Error found at double precision =        Poly: 1.441088e-16   Cheb: 1.864792e-16


      static const T Y = 1.03937530517578125F;
      static const T P[9] = {    
         -0.00356165148914447597995L,
         -0.0229930320357982333406L,
         -0.0449814350482277917716L,
         -0.0453759383048193402336L,
         -0.0272050837209380717069L,
         -0.00994403059883350813295L,
         -0.00207592267812291726961L,
         -0.000192178045857733706044L,
         -0.113161784705911400295e-9L
      };
      static const T Q[9] = {    
         1L,
         2.84354408840148561131L,
         3.6599610090072393012L,
         2.75088464344293083595L,
         1.2985244073998398643L,
         0.383213198510794507409L,
         0.0651165455496281337831L,
         0.00488071077519227853585L
      };
      T t = z / 10 - 3;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      result *= exp(z) / z;
      result += z;
   }
   else
   {
      // Max Error found at double precision =        3.381886e-17
      static const T exp40 = static_cast<T>(2.35385266837019985407899910749034804508871617254555467236651e17L);
      static const T Y= 1.013065338134765625F;
      static const T P[6] = {    
         -0.0130653381347656243849L,
         0.19029710559486576682L,
         94.7365094537197236011L,
         -2516.35323679844256203L,
         18932.0850014925993025L,
         -38703.1431362056714134L
      };
      static const T Q[7] = {    
         1L,
         61.9733592849439884145L,
         -2354.56211323420194283L,
         22329.1459489893079041L,
         -70126.245140396567133L,
         54738.2833147775537106L,
         8297.16296356518409347L
      };
      T t = 1 / z;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      if(z < 41)
         result *= exp(z) / z;
      else
      {
         // Avoid premature overflow if we can:
         t = z - 40;
         if(t > tools::log_max_value<T>())
         {
            result = policies::raise_overflow_error<T>(function, 0, pol);
         }
         else
         {
            result *= exp(z - 40) / z;
            if(result > tools::max_value<T>() / exp40)
            {
               result = policies::raise_overflow_error<T>(function, 0, pol);
            }
            else
            {
               result *= exp40;
            }
         }
      }
      result += z;
   }
   return result;
}

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

   T result;

   if(z <= 6)
   {
      // Maximum Deviation Found:                     3.883e-21
      // Expected Error Term:                         3.883e-21
      // Max Error found at long double precision =   Poly: 3.344801e-19   Cheb: 4.989937e-19

      static const T P[11] = {    
         2.98677224343598593764L,
         0.25891613550886736592L,
         0.789323584998672832285L,
         0.092432587824602399339L,
         0.0514236978728625906656L,
         0.00658477469745132977921L,
         0.00124914538197086254233L,
         0.000131429679565472408551L,
         0.11293331317982763165e-4L,
         0.629499283139417444244e-6L,
         0.177833045143692498221e-7L
      };
      static const T Q[9] = {    
         1L,