expint.hpp

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

//  Copyright John Maddock 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_EXPINT_HPP
#define BOOST_MATH_EXPINT_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/math/tools/fraction.hpp>
#include <boost/math/tools/series.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/pow.hpp>

namespace boost{ namespace math{

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   expint(unsigned n, T z, const Policy& /*pol*/);
   
namespace detail{

template <class T>
inline T expint_1_rational(const T& z, const mpl::int_<0>&)
{
   // this function is never actually called
   BOOST_ASSERT(0);
   return z;
}

template <class T>
T expint_1_rational(const T& z, const mpl::int_<53>&)
{
   BOOST_MATH_STD_USING
   T result;
   if(z <= 1)
   {
      // Maximum Deviation Found:                     2.006e-18
      // Expected Error Term:                         2.006e-18
      // Max error found at double precision:         2.760e-17
      static const T Y = 0.66373538970947265625F;
      static const T P[6] = {    
         0.0865197248079397976498L,
         0.0320913665303559189999L,
         -0.245088216639761496153L,
         -0.0368031736257943745142L,
         -0.00399167106081113256961L,
         -0.000111507792921197858394L
      };
      static const T Q[6] = {    
         1L,
         0.37091387659397013215L,
         0.056770677104207528384L,
         0.00427347600017103698101L,
         0.000131049900798434683324L,
         -0.528611029520217142048e-6L
      };
      result = tools::evaluate_polynomial(P, z) 
         / tools::evaluate_polynomial(Q, z);
      result += z - log(z) - Y;
   }
   else if(z < -boost::math::tools::log_min_value<T>())
   {
      // Maximum Deviation Found (interpolated):      1.444e-17
      // Max error found at double precision:         3.119e-17
      static const T P[11] = {    
         -0.121013190657725568138e-18L,
         -0.999999999999998811143L,
         -43.3058660811817946037L,
         -724.581482791462469795L,
         -6046.8250112711035463L,
         -27182.6254466733970467L,
         -66598.2652345418633509L,
         -86273.1567711649528784L,
         -54844.4587226402067411L,
         -14751.4895786128450662L,
         -1185.45720315201027667L
      };
      static const T Q[12] = {    
         1L,
         45.3058660811801465927L,
         809.193214954550328455L,
         7417.37624454689546708L,
         38129.5594484818471461L,
         113057.05869159631492L,
         192104.047790227984431L,
         180329.498380501819718L,
         86722.3403467334749201L,
         18455.4124737722049515L,
         1229.20784182403048905L,
         -0.776491285282330997549L
      };
      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>
T expint_1_rational(const T& z, const mpl::int_<64>&)
{
   BOOST_MATH_STD_USING
   T result;
   if(z <= 1)
   {
      // Maximum Deviation Found:                     3.807e-20
      // Expected Error Term:                         3.807e-20
      // Max error found at long double precision:    6.249e-20

      static const T Y = 0.66373538970947265625F;
      static const T P[6] = {    
         0.0865197248079397956816L,
         0.0275114007037026844633L,
         -0.246594388074877139824L,
         -0.0237624819878732642231L,
         -0.00259113319641673986276L,
         0.30853660894346057053e-4L
      };
      static const T Q[7] = {    
         1L,
         0.317978365797784100273L,
         0.0393622602554758722511L,
         0.00204062029115966323229L,
         0.732512107100088047854e-5L,
         -0.202872781770207871975e-5L,
         0.52779248094603709945e-7L
      };
      result = tools::evaluate_polynomial(P, z) 
         / tools::evaluate_polynomial(Q, z);
      result += z - log(z) - Y;
   }
   else if(z < -boost::math::tools::log_min_value<T>())
   {
      // Maximum Deviation Found (interpolated):     2.220e-20
      // Max error found at long double precision:   1.346e-19
      static const T P[14] = {    
         -0.534401189080684443046e-23L,
         -0.999999999999999999905L,
         -62.1517806091379402505L,
         -1568.45688271895145277L,
         -21015.3431990874009619L,
         -164333.011755931661949L,
         -777917.270775426696103L,
         -2244188.56195255112937L,
         -3888702.98145335643429L,
         -3909822.65621952648353L,
         -2149033.9538897398457L,
         -584705.537139793925189L,
         -65815.2605361889477244L,
         -2038.82870680427258038L
      };
      static const T Q[14] = {    
         1L,
         64.1517806091379399478L,
         1690.76044393722763785L,
         24035.9534033068949426L,
         203679.998633572361706L,
         1074661.58459976978285L,
         3586552.65020899358773L,
         7552186.84989547621411L,
         9853333.79353054111434L,
         7689642.74550683631258L,
         3385553.35146759180739L,
         763218.072732396428725L,
         73930.2995984054930821L,
         2063.86994219629165937L
      };
      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>
T expint_1_rational(const T& z, const mpl::int_<113>&)
{
   BOOST_MATH_STD_USING
   T result;
   if(z <= 1)
   {
      // Maximum Deviation Found:                     2.477e-35
      // Expected Error Term:                         2.477e-35
      // Max error found at long double precision:    6.810e-35

      static const T Y = 0.66373538970947265625F;
      static const T P[10] = {    
         0.0865197248079397956434879099175975937L,
         0.0369066175910795772830865304506087759L,
         -0.24272036838415474665971599314725545L,
         -0.0502166331248948515282379137550178307L,
         -0.00768384138547489410285101483730424919L,
         -0.000612574337702109683505224915484717162L,
         -0.380207107950635046971492617061708534e-4L,
         -0.136528159460768830763009294683628406e-5L,
         -0.346839106212658259681029388908658618e-7L,
         -0.340500302777838063940402160594523429e-9L
      };
      static const T Q[10] = {    
         1L,
         0.426568827778942588160423015589537302L,
         0.0841384046470893490592450881447510148L,
         0.0100557215850668029618957359471132995L,
         0.000799334870474627021737357294799839363L,
         0.434452090903862735242423068552687688e-4L,
         0.15829674748799079874182885081231252e-5L,
         0.354406206738023762100882270033082198e-7L,
         0.369373328141051577845488477377890236e-9L,
         -0.274149801370933606409282434677600112e-12L
      };
      result = tools::evaluate_polynomial(P, z) 
         / tools::evaluate_polynomial(Q, z);
      result += z - log(z) - Y;
   }
   else if(z <= 4)
   {
      // Max error in interpolated form:             5.614e-35
      // Max error found at long double precision:   7.979e-35

      static const T Y = 0.70190334320068359375F;

      static const T P[17] = {    
         0.298096656795020369955077350585959794L,
         12.9314045995266142913135497455971247L,
         226.144334921582637462526628217345501L,
         2070.83670924261732722117682067381405L,
         10715.1115684330959908244769731347186L,
         30728.7876355542048019664777316053311L,
         38520.6078609349855436936232610875297L,
         -27606.0780981527583168728339620565165L,
         -169026.485055785605958655247592604835L,
         -254361.919204983608659069868035092282L,
         -195765.706874132267953259272028679935L,
         -83352.6826013533205474990119962408675L,
         -19251.6828496869586415162597993050194L,
         -2226.64251774578542836725386936102339L,
         -109.009437301400845902228611986479816L,
         -1.51492042209561411434644938098833499L
      };
      static const T Q[16] = {    
         1L,
         46.734521442032505570517810766704587L,
         908.694714348462269000247450058595655L,
         9701.76053033673927362784882748513195L,
         63254.2815292641314236625196594947774L,
         265115.641285880437335106541757711092L,
         732707.841188071900498536533086567735L,
         1348514.02492635723327306628712057794L,
         1649986.81455283047769673308781585991L,
         1326000.828522976970116271208812099L,
         683643.09490612171772350481773951341L,
         217640.505137263607952365685653352229L,
         40288.3467237411710881822569476155485L,
         3932.89353979531632559232883283175754L,
         169.845369689596739824177412096477219L,
         2.17607292280092201170768401876895354L
      };
      T recip = 1 / z;
      result = Y + tools::evaluate_polynomial(P, recip)
         / tools::evaluate_polynomial(Q, recip);
      result *= exp(-z) * recip;
   }
   else if(z < -boost::math::tools::log_min_value<T>())
   {
      // Max error in interpolated form:             4.413e-35
      // Max error found at long double precision:   8.928e-35

      static const T P[19] = {    
         -0.559148411832951463689610809550083986e-40L,
         -0.999999999999999999999999999999999997L,
         -166.542326331163836642960118190147367L,
         -12204.639128796330005065904675153652L,
         -520807.069767086071806275022036146855L,
         -14435981.5242137970691490903863125326L,
         -274574945.737064301247496460758654196L,
         -3691611582.99810039356254671781473079L,
         -35622515944.8255047299363690814678763L,
         -248040014774.502043161750715548451142L,
         -1243190389769.53458416330946622607913L,
         -4441730126135.54739052731990368425339L,
         -11117043181899.7388524310281751971366L,
         -18976497615396.9717776601813519498961L,
         -21237496819711.1011661104761906067131L,
         -14695899122092.5161620333466757812848L,
         -5737221535080.30569711574295785864903L,
         -1077042281708.42654526404581272546244L,
         -68028222642.1941480871395695677675137L
      };
      static const T Q[20] = {    
         1L,
         168.542326331163836642960118190147311L,
         12535.7237814586576783518249115343619L,
         544891.263372016404143120911148640627L,
         15454474.7241010258634446523045237762L,
         302495899.896629522673410325891717381L,
         4215565948.38886507646911672693270307L,
         42552409471.7951815668506556705733344L,
         313592377066.753173979584098301610186L,
         1688763640223.4541980740597514904542L,
         6610992294901.59589748057620192145704L,
         18601637235659.6059890851321772682606L,
         36944278231087.2571020964163402941583L,
         50425858518481.7497071917028793820058L,
         45508060902865.0899967797848815980644L,
         25649955002765.3817331501988304758142L,
         8259575619094.6518520988612711292331L,
         1299981487496.12607474362723586264515L,
         70242279152.8241187845178443118302693L,
         -37633302.9409263839042721539363416685L
      };
      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, 
      tools::digits<T>(),
      max_iter);
   policies::check_series_iterations("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