expint.hpp

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

         -1.20352377969742325748L,
         0.66707904942606479811L,
         -0.223014531629140771914L,
         0.0493340022262908008636L,
         -0.00741934273050807310677L,
         0.00074353567782087939294L,
         -0.455861727069603367656e-4L,
         0.131515429329812837701e-5L
      };

      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:                     2.622e-21
      // Expected Error Term:                         -2.622e-21
      // Max Error found at long double precision =   Poly: 1.208328e-20   Cheb: 1.073723e-20

      static const T Y = 1.158985137939453125F;
      static const T P[9] = {    
         0.00139324086199409049399L,
         -0.0345238388952337563247L,
         -0.0382065278072592940767L,
         -0.0156117003070560727392L,
         -0.00383276012430495387102L,
         -0.000697070540945496497992L,
         -0.877310384591205930343e-4L,
         -0.623067256376494930067e-5L,
         -0.377246883283337141444e-6L
      };
      static const T Q[10] = {    
         1L,
         1.08073635708902053767L,
         0.553681133533942532909L,
         0.176763647137553797451L,
         0.0387891748253869928121L,
         0.0060603004848394727017L,
         0.000670519492939992806051L,
         0.4947357050100855646e-4L,
         0.204339282037446434827e-5L,
         0.146951181174930425744e-7L
      };
      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:                     3.220e-20
      // Expected Error Term:                         3.220e-20
      // Max Error found at long double precision =   Poly: 7.696841e-20   Cheb: 6.205163e-20


      static const T Y = 1.0869731903076171875F;
      static const T P[10] = {    
         -0.00893891094356946995368L,
         -0.0487562980088748775943L,
         -0.0670568657950041926085L,
         -0.0509577352851442932713L,
         -0.02551800927409034206L,
         -0.00892913759760086687083L,
         -0.00224469630207344379888L,
         -0.000392477245911296982776L,
         -0.44424044184395578775e-4L,
         -0.252788029251437017959e-5L
      };
      static const T Q[10] = {    
         1L,
         2.00323265503572414261L,
         1.94688958187256383178L,
         1.19733638134417472296L,
         0.513137726038353385661L,
         0.159135395578007264547L,
         0.0358233587351620919881L,
         0.0056716655597009417875L,
         0.000577048986213535829925L,
         0.290976943033493216793e-4L
      };
      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:                     2.940e-21
      // Expected Error Term:                         -2.938e-21
      // Max Error found at long double precision =   Poly: 3.419893e-19   Cheb: 3.359874e-19

      static const T Y = 1.03937530517578125F;
      static const T P[12] = {    
         -0.00356165148914447278177L,
         -0.0240235006148610849678L,
         -0.0516699967278057976119L,
         -0.0586603078706856245674L,
         -0.0409960120868776180825L,
         -0.0185485073689590665153L,
         -0.00537842101034123222417L,
         -0.000920988084778273760609L,
         -0.716742618812210980263e-4L,
         -0.504623302166487346677e-9L,
         0.712662196671896837736e-10L,
         -0.533769629702262072175e-11L
      };
      static const T Q[9] = {    
         1L,
         3.13286733695729715455L,
         4.49281223045653491929L,
         3.84900294427622911374L,
         2.15205199043580378211L,
         0.802912186540269232424L,
         0.194793170017818925388L,
         0.0280128013584653182994L,
         0.00182034930799902922549L
      };
      T t = z / 10 - 3;
      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
   {
      // Maximum Deviation Found:                     3.536e-20
      // Max Error found at long double precision =   Poly: 1.310671e-19   Cheb: 8.630943e-11

      static const T exp40 = static_cast<T>(2.35385266837019985407899910749034804508871617254555467236651e17L);
      static const T Y= 1.013065338134765625F;
      static const T P[9] = {    
         -0.0130653381347656250004L,
         0.644487780349757303739L,
         143.995670348227433964L,
         -13918.9322758014173709L,
         476260.975133624194484L,
         -7437102.15135982802122L,
         53732298.8764767916542L,
         -160695051.957997452509L,
         137839271.592778020028L
      };
      static const T Q[9] = {    
         1L,
         27.2103343964943718802L,
         -8785.48528692879413676L,
         397530.290000322626766L,
         -7356441.34957799368252L,
         63050914.5343400957524L,
         -246143779.638307701369L,
         384647824.678554961174L,
         -166288297.874583961493L
      };
      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_<113>& 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:                     1.230e-36
      // Expected Error Term:                         -1.230e-36
      // Max Error found at long double precision =   Poly: 4.355299e-34   Cheb: 7.512581e-34


      static const T P[15] = {    
         2.98677224343598593765287235997328555L,
         -0.333256034674702967028780537349334037L,
         0.851831522798101228384971644036708463L,
         -0.0657854833494646206186773614110374948L,
         0.0630065662557284456000060708977935073L,
         -0.00311759191425309373327784154659649232L,
         0.00176213568201493949664478471656026771L,
         -0.491548660404172089488535218163952295e-4L,
         0.207764227621061706075562107748176592e-4L,
         -0.225445398156913584846374273379402765e-6L,
         0.996939977231410319761273881672601592e-7L,
         0.212546902052178643330520878928100847e-9L,
         0.154646053060262871360159325115980023e-9L,
         0.143971277122049197323415503594302307e-11L,
         0.306243138978114692252817805327426657e-13L
      };
      static const T Q[15] = {    
         1L,
         -1.40178870313943798705491944989231793L,
         0.943810968269701047641218856758605284L,
         -0.405026631534345064600850391026113165L,
         0.123924153524614086482627660399122762L,
         -0.0286364505373369439591132549624317707L,
         0.00516148845910606985396596845494015963L,
         -0.000738330799456364820380739850924783649L,
         0.843737760991856114061953265870882637e-4L,
         -0.767957673431982543213661388914587589e-5L,
         0.549136847313854595809952100614840031e-6L,
         -0.299801381513743676764008325949325404e-7L,
         0.118419479055346106118129130945423483e-8L,
         -0.30372295663095470359211949045344607e-10L,
         0.382742953753485333207877784720070523e-12L
      };

      static const T r1 = static_cast<T>(1677624236387711.0L) / 4503599627370496uLL;
      static const T r2 = (static_cast<T>(266514582277687.0L) / 4503599627370496uLL) / 4503599627370496uLL;
      static const T r3 = static_cast<T>(0.283806480836357377069325311780969887585024578164571984232357e-31L);
      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) - r3;
      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:                     7.779e-36
      // Expected Error Term:                         -7.779e-36
      // Max Error found at long double precision =   Poly: 2.576723e-35   Cheb: 1.236001e-34

      static const T Y = 1.158985137939453125F;
      static const T P[15] = {    
         0.00139324086199409049282472239613554817L,
         -0.0338173111691991289178779840307998955L,
         -0.0555972290794371306259684845277620556L,
         -0.0378677976003456171563136909186202177L,
         -0.0152221583517528358782902783914356667L,
         -0.00428283334203873035104248217403126905L,
         -0.000922782631491644846511553601323435286L,
         -0.000155513428088853161562660696055496696L,
         -0.205756580255359882813545261519317096e-4L,
         -0.220327406578552089820753181821115181e-5L,
         -0.189483157545587592043421445645377439e-6L,
         -0.122426571518570587750898968123803867e-7L,
         -0.635187358949437991465353268374523944e-9L,
         -0.203015132965870311935118337194860863e-10L,
         -0.384276705503357655108096065452950822e-12L
      };
      static const T Q[15] = {    
         1L,
         1.58784732785354597996617046880946257L,
         1.18550755302279446339364262338114098L,
         0.55598993549661368604527040349702836L,
         0.184290888380564236919107835030984453L,
         0.0459658051803613282360464632326866113L,
         0.0089505064268613225167835599456014705L,
         0.00139042673882987693424772855926289077L,
         0.000174210708041584097450805790176479012L,
         0.176324034009707558089086875136647376e-4L,
         0.142935845999505649273084545313710581e-5L,
         0.907502324487057260675816233312747784e-7L,
         0.431044337808893270797934621235918418e-8L,
         0.139007266881450521776529705677086902e-9L,
         0.234715286125516430792452741830364672e-11L
      };
      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 <= 18)
   {
      // Maximum Deviation Found:                     1.082e-34
      // Expected Error Term:                         1.080e-34
      // Max Error found at long double precision =   Poly: 1.958294e-34   Cheb: 2.472261e-34


      static const T Y = 1.091579437255859375F;
      static const T P[17] = {    
         -0.00685089599550151282724924894258520532L,
         -0.0443313550253580053324487059748497467L,
         -0.071538561252424027443296958795814874L,
         -0.0622923153354102682285444067843300583L,
         -0.0361631270264607478205393775461208794L,
         -0.0153192826839624850298106509601033261L,
         -0.00496967904961260031539602977748408242L,
         -0.00126989079663425780800919171538920589L,
         -0.000258933143097125199914724875206326698L,
         -0.422110326689204794443002330541441956e-4L,
         -0.546004547590412661451073996127115221e-5L,
         -0.546775260262202177131068692199272241e-6L,
         -0.404157632825805803833379568956559215e-7L,
         -0.200612596196561323832327013027419284e-8L,
         -0.502538501472133913417609379765434153e-10L,
         -0.326283053716799774936661568391296584e-13L,
         0.869226483473172853557775877908693647e-15L
      };
      static const T Q[15] = {    
         1L,
         2.23227220874479061894038229141871087L,
         2.40221000361027971895657505660959863L,
         1.65476320985936174728238416007084214L,
         0.816828602963895720369875535001248227L,
         0.306337922909446903672123418670921066L,
         0.0902400121654409267774593230720600752L,
         0.0212708882169429206498765100993228086L,
         0.00404442626252467471957713495828165491L,
         0.0006195601618842253612635241404054589L,
         0.755930932686543009521454653994321843e-4L,
         0.716004532773778954193609582677482803e-5L,
         0.500881663076471627699290821742924233e-6L,
         0.233593219218823384508105943657387644e-7L,
         0.554900353169148897444104962034267682e-9L
      };
      T t = z / 4 - 3.5;
      result = Y + tools::evaluate_polynomial(P, t)
         / tools::evaluate_polynomial(Q, t);
      result *= exp(z) / z;
      result += z;
   }
   else if(z <= 26)
   {
      // Maximum Deviation Found:                     3.163e-35
      // Expected Error Term:                         3.163e-35
      // Max Error found at long double precision =   Poly: 4.158110e-35   Cheb: 5.385532e-35

      static const T Y = 1.051731109619140625F;
      static const T P[14] = {