expint.hpp

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

         -0.00144552494420652573815404828020593565L,
         -0.0126747451594545338365684731262912741L,
         -0.01757394877502366717526779263438073L,
         -0.0126838952395506921945756139424722588L,
         -0.0060045057928894974954756789352443522L,
         -0.00205349237147226126653803455793107903L,
         -0.000532606040579654887676082220195624207L,
         -0.000107344687098019891474772069139014662L,
         -0.169536802705805811859089949943435152e-4L,
         -0.20863311729206543881826553010120078e-5L,
         -0.195670358542116256713560296776654385e-6L,
         -0.133291168587253145439184028259772437e-7L,
         -0.595500337089495614285777067722823397e-9L,
         -0.133141358866324100955927979606981328e-10L
      };
      static const T Q[14] = {    
         1L,
         1.72490783907582654629537013560044682L,
         1.44524329516800613088375685659759765L,
         0.778241785539308257585068744978050181L,
         0.300520486589206605184097270225725584L,
         0.0879346899691339661394537806057953957L,
         0.0200802415843802892793583043470125006L,
         0.00362842049172586254520256100538273214L,
         0.000519731362862955132062751246769469957L,
         0.584092147914050999895178697392282665e-4L,
         0.501851497707855358002773398333542337e-5L,
         0.313085677467921096644895738538865537e-6L,
         0.127552010539733113371132321521204458e-7L,
         0.25737310826983451144405899970774587e-9L
      };
      T t = z / 4 - 5.5;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result *= exp(z) / z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result += z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
   }
   else if(z <= 42)
   {
      // Maximum Deviation Found:                     7.972e-36
      // Expected Error Term:                         7.962e-36
      // Max Error found at long double precision =   Poly: 1.711721e-34   Cheb: 3.100018e-34

      static const T Y = 1.032726287841796875F;
      static const T P[15] = {    
         -0.00141056919297307534690895009969373233L,
         -0.0123384175302540291339020257071411437L,
         -0.0298127270706864057791526083667396115L,
         -0.0390686759471630584626293670260768098L,
         -0.0338226792912607409822059922949035589L,
         -0.0211659736179834946452561197559654582L,
         -0.0100428887460879377373158821400070313L,
         -0.00370717396015165148484022792801682932L,
         -0.0010768667551001624764329000496561659L,
         -0.000246127328761027039347584096573123531L,
         -0.437318110527818613580613051861991198e-4L,
         -0.587532682329299591501065482317771497e-5L,
         -0.565697065670893984610852937110819467e-6L,
         -0.350233957364028523971768887437839573e-7L,
         -0.105428907085424234504608142258423505e-8L
      };
      static const T Q[16] = {    
         1L,
         3.17261315255467581204685605414005525L,
         4.85267952971640525245338392887217426L,
         4.74341914912439861451492872946725151L,
         3.31108463283559911602405970817931801L,
         1.74657006336994649386607925179848899L,
         0.718255607416072737965933040353653244L,
         0.234037553177354542791975767960643864L,
         0.0607470145906491602476833515412605389L,
         0.0125048143774226921434854172947548724L,
         0.00201034366420433762935768458656609163L,
         0.000244823338417452367656368849303165721L,
         0.213511655166983177960471085462540807e-4L,
         0.119323998465870686327170541547982932e-5L,
         0.322153582559488797803027773591727565e-7L,
         -0.161635525318683508633792845159942312e-16L
      };
      T t = z / 8 - 4.25;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result *= exp(z) / z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result += z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
   }
   else if(z <= 56)
   {
      // Maximum Deviation Found:                     4.469e-36
      // Expected Error Term:                         4.468e-36
      // Max Error found at long double precision =   Poly: 1.288958e-35   Cheb: 2.304586e-35

      static const T Y = 1.0216197967529296875F;
      static const T P[12] = {    
         -0.000322999116096627043476023926572650045L,
         -0.00385606067447365187909164609294113346L,
         -0.00686514524727568176735949971985244415L,
         -0.00606260649593050194602676772589601799L,
         -0.00334382362017147544335054575436194357L,
         -0.00126108534260253075708625583630318043L,
         -0.000337881489347846058951220431209276776L,
         -0.648480902304640018785370650254018022e-4L,
         -0.87652644082970492211455290209092766e-5L,
         -0.794712243338068631557849449519994144e-6L,
         -0.434084023639508143975983454830954835e-7L,
         -0.107839681938752337160494412638656696e-8L
      };
      static const T Q[12] = {    
         1L,
         2.09913805456661084097134805151524958L,
         2.07041755535439919593503171320431849L,
         1.26406517226052371320416108604874734L,
         0.529689923703770353961553223973435569L,
         0.159578150879536711042269658656115746L,
         0.0351720877642000691155202082629857131L,
         0.00565313621289648752407123620997063122L,
         0.000646920278540515480093843570291218295L,
         0.499904084850091676776993523323213591e-4L,
         0.233740058688179614344680531486267142e-5L,
         0.498800627828842754845418576305379469e-7L
      };
      T t = z / 7 - 7;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result *= exp(z) / z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result += z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
   }
   else if(z <= 84)
   {
      // Maximum Deviation Found:                     5.588e-35
      // Expected Error Term:                         -5.566e-35
      // Max Error found at long double precision =   Poly: 9.976345e-35   Cheb: 8.358865e-35

      static const T Y = 1.015148162841796875F;
      static const T P[11] = {    
         -0.000435714784725086961464589957142615216L,
         -0.00432114324353830636009453048419094314L,
         -0.0100740363285526177522819204820582424L,
         -0.0116744115827059174392383504427640362L,
         -0.00816145387784261141360062395898644652L,
         -0.00371380272673500791322744465394211508L,
         -0.00112958263488611536502153195005736563L,
         -0.000228316462389404645183269923754256664L,
         -0.29462181955852860250359064291292577e-4L,
         -0.21972450610957417963227028788460299e-5L,
         -0.720558173805289167524715527536874694e-7L
      };
      static const T Q[11] = {    
         1L,
         2.95918362458402597039366979529287095L,
         3.96472247520659077944638411856748924L,
         3.15563251550528513747923714884142131L,
         1.64674612007093983894215359287448334L,
         0.58695020129846594405856226787156424L,
         0.144358385319329396231755457772362793L,
         0.024146911506411684815134916238348063L,
         0.0026257132337460784266874572001650153L,
         0.000167479843750859222348869769094711093L,
         0.475673638665358075556452220192497036e-5L
      };
      T t = z / 14 - 5;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result *= exp(z) / z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
      result += z;
      BOOST_MATH_INSTRUMENT_VARIABLE(result)
   }
   else if(z <= 210)
   {
      // Maximum Deviation Found:                     4.448e-36
      // Expected Error Term:                         4.445e-36
      // Max Error found at long double precision =   Poly: 2.058532e-35   Cheb: 2.165465e-27

      static const T Y= 1.00849151611328125F;
      static const T P[9] = {    
         -0.0084915161132812500000001440233607358L,
         1.84479378737716028341394223076147872L,
         -130.431146923726715674081563022115568L,
         4336.26945491571504885214176203512015L,
         -76279.0031974974730095170437591004177L,
         729577.956271997673695191455111727774L,
         -3661928.69330208734947103004900349266L,
         8570600.041606912735872059184527855L,
         -6758379.93672362080947905580906028645L
      };
      static const T Q[10] = {    
         1L,
         -99.4868026047611434569541483506091713L,
         3879.67753690517114249705089803055473L,
         -76495.82413252517165830203774900806L,
         820773.726408311894342553758526282667L,
         -4803087.64956923577571031564909646579L,
         14521246.227703545012713173740895477L,
         -19762752.0196769712258527849159393044L,
         8354144.67882768405803322344185185517L,
         355076.853106511136734454134915432571L
      };
      T t = 1 / z;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      result *= exp(z) / z;
      result += z;
   }
   else // z > 210
   {
      // Maximum Deviation Found:                     3.963e-37
      // Expected Error Term:                         3.963e-37
      // Max Error found at long double precision =   Poly: 1.248049e-36   Cheb: 2.843486e-29

      static const T exp40 = static_cast<T>(2.35385266837019985407899910749034804508871617254555467236651e17L);
      static const T Y= 1.00252532958984375F;
      static const T P[8] = {    
         -0.00252532958984375000000000000000000085L,
         1.16591386866059087390621952073890359L,
         -67.8483431314018462417456828499277579L,
         1567.68688154683822956359536287575892L,
         -17335.4683325819116482498725687644986L,
         93632.6567462673524739954389166550069L,
         -225025.189335919133214440347510936787L,
         175864.614717440010942804684741336853L
      };
      static const T Q[9] = {    
         1L,
         -65.6998869881600212224652719706425129L,
         1642.73850032324014781607859416890077L,
         -19937.2610222467322481947237312818575L,
         124136.267326632742667972126625064538L,
         -384614.251466704550678760562965502293L,
         523355.035910385688578278384032026998L,
         -217809.552260834025885677791936351294L,
         -8555.81719551123640677261226549550872L
      };
      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>
inline typename tools::promote_args<T>::type
   expint_forwarder(T z, const Policy& /*pol*/, mpl::true_ const&)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   typedef typename policies::precision<result_type, Policy>::type precision_type;
   typedef typename policies::normalise<
      Policy, 
      policies::promote_float<false>, 
      policies::promote_double<false>, 
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;
   typedef typename mpl::if_<
      mpl::less_equal<precision_type, mpl::int_<0> >,
      mpl::int_<0>,
      typename mpl::if_<
         mpl::less_equal<precision_type, mpl::int_<53> >,
         mpl::int_<53>,  // double
         typename mpl::if_<
            mpl::less_equal<precision_type, mpl::int_<64> >,
            mpl::int_<64>, // 80-bit long double
            typename mpl::if_<
               mpl::less_equal<precision_type, mpl::int_<113> >,
               mpl::int_<113>, // 128-bit long double
               mpl::int_<0> // too many bits, use generic version.
            >::type
         >::type
      >::type
   >::type tag_type;

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_i_imp(
      static_cast<value_type>(z),
      forwarding_policy(),
      tag_type()), "boost::math::expint<%1%>(%1%)");
}

template <class T>
inline typename tools::promote_args<T>::type
expint_forwarder(unsigned n, T z, const mpl::false_&)
{
   return boost::math::expint(n, z, policies::policy<>());
}

} // namespace detail

template <class T, class Policy>
inline typename tools::promote_args<T>::type
   expint(unsigned n, T z, const Policy& /*pol*/)
{
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   typedef typename policies::precision<result_type, Policy>::type precision_type;
   typedef typename policies::normalise<
      Policy, 
      policies::promote_float<false>, 
      policies::promote_double<false>, 
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;
   typedef typename mpl::if_<
      mpl::less_equal<precision_type, mpl::int_<0> >,
      mpl::int_<0>,
      typename mpl::if_<
         mpl::less_equal<precision_type, mpl::int_<53> >,
         mpl::int_<53>,  // double
         typename mpl::if_<
            mpl::less_equal<precision_type, mpl::int_<64> >,
            mpl::int_<64>, // 80-bit long double
            typename mpl::if_<
               mpl::less_equal<precision_type, mpl::int_<113> >,
               mpl::int_<113>, // 128-bit long double
               mpl::int_<0> // too many bits, use generic version.
            >::type
         >::type
      >::type
   >::type tag_type;

   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_imp(
      n,
      static_cast<value_type>(z),
      forwarding_policy(),
      tag_type()), "boost::math::expint<%1%>(unsigned, %1%)");
}

template <class T, class U>
inline typename detail::expint_result<T, U>::type
   expint(T const z, U const u)
{
   typedef typename policies::is_policy<U>::type tag_type;
   return detail::expint_forwarder(z, u, tag_type());
}

template <class T>
inline typename tools::promote_args<T>::type
   expint(T z)
{
   return expint(z, policies::policy<>());
}

}} // namespaces

#endif // BOOST_MATH_EXPINT_HPP