erf.hpp

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

            0.00656970902163248872837262539337601845L,
            0.00120282643299089441390490459256235021L,
            0.000142123229065182650020762792081622986L,
            0.991531438367015135346716277792989347e-5L,
            0.312857043762117596999398067153076051e-6L,
         };
         static const T Q[] = {    
            1L,
            2.13506082409097783827103424943508554L,
            2.06399257267556230937723190496806215L,
            1.18678481279932541314830499880691109L,
            0.447733186643051752513538142316799562L,
            0.11505680005657879437196953047542148L,
            0.020163993632192726170219663831914034L,
            0.00232708971840141388847728782209730585L,
            0.000160733201627963528519726484608224112L,
            0.507158721790721802724402992033269266e-5L,
            0.18647774409821470950544212696270639e-12L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
         result *= exp(-z * z) / z;
      }
      else if (z < 3)
      {
         // Maximum Deviation Found:                     3.575e-36
         // Expected Error Term:                         3.575e-36
         // Maximum Relative Change in Control Points:   7.103e-05
         // Max Error found at long double precision =   5.794737e-36
         static const T Y = 0.52896785736083984375f;
         static const T P[] = {    
            -0.00902152521745813634562524098263360074L,
            0.0145207142776691539346923710537580927L,
            0.0301681239582193983824211995978678571L,
            0.0215548540823305814379020678660434461L,
            0.00864683476267958365678294164340749949L,
            0.00219693096885585491739823283511049902L,
            0.000364961639163319762492184502159894371L,
            0.388174251026723752769264051548703059e-4L,
            0.241918026931789436000532513553594321e-5L,
            0.676586625472423508158937481943649258e-7L,
         };
         static const T Q[] = {    
            1L,
            1.93669171363907292305550231764920001L,
            1.69468476144051356810672506101377494L,
            0.880023580986436640372794392579985511L,
            0.299099106711315090710836273697708402L,
            0.0690593962363545715997445583603382337L,
            0.0108427016361318921960863149875360222L,
            0.00111747247208044534520499324234317695L,
            0.686843205749767250666787987163701209e-4L,
            0.192093541425429248675532015101904262e-5L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 2.25f) / tools::evaluate_polynomial(Q, z - 2.25f);
         result *= exp(-z * z) / z;
      }
      else if(z < 3.5)
      {
         // Maximum Deviation Found:                     8.126e-37
         // Expected Error Term:                         -8.126e-37
         // Maximum Relative Change in Control Points:   1.363e-04
         // Max Error found at long double precision =   1.747062e-36
         static const T Y = 0.54037380218505859375f;
         static const T P[] = {    
            -0.0033703486408887424921155540591370375L,
            0.0104948043110005245215286678898115811L,
            0.0148530118504000311502310457390417795L,
            0.00816693029245443090102738825536188916L,
            0.00249716579989140882491939681805594585L,
            0.0004655591010047353023978045800916647L,
            0.531129557920045295895085236636025323e-4L,
            0.343526765122727069515775194111741049e-5L,
            0.971120407556888763695313774578711839e-7L,
         };
         static const T Q[] = {    
            1L,
            1.59911256167540354915906501335919317L,
            1.136006830764025173864831382946934L,
            0.468565867990030871678574840738423023L,
            0.122821824954470343413956476900662236L,
            0.0209670914950115943338996513330141633L,
            0.00227845718243186165620199012883547257L,
            0.000144243326443913171313947613547085553L,
            0.407763415954267700941230249989140046e-5L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 3.0f) / tools::evaluate_polynomial(Q, z - 3.0f);
         result *= exp(-z * z) / z;
      }
      else if(z < 5.5)
      {
         // Maximum Deviation Found:                     5.804e-36
         // Expected Error Term:                         -5.803e-36
         // Maximum Relative Change in Control Points:   2.475e-05
         // Max Error found at long double precision =   1.349545e-35
         static const T Y = 0.55000019073486328125f;
         static const T P[] = {    
            0.00118142849742309772151454518093813615L,
            0.0072201822885703318172366893469382745L,
            0.0078782276276860110721875733778481505L,
            0.00418229166204362376187593976656261146L,
            0.00134198400587769200074194304298642705L,
            0.000283210387078004063264777611497435572L,
            0.405687064094911866569295610914844928e-4L,
            0.39348283801568113807887364414008292e-5L,
            0.248798540917787001526976889284624449e-6L,
            0.929502490223452372919607105387474751e-8L,
            0.156161469668275442569286723236274457e-9L,
         };
         static const T Q[] = {    
            1L,
            1.52955245103668419479878456656709381L,
            1.06263944820093830054635017117417064L,
            0.441684612681607364321013134378316463L,
            0.121665258426166960049773715928906382L,
            0.0232134512374747691424978642874321434L,
            0.00310778180686296328582860464875562636L,
            0.000288361770756174705123674838640161693L,
            0.177529187194133944622193191942300132e-4L,
            0.655068544833064069223029299070876623e-6L,
            0.11005507545746069573608988651927452e-7L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 4.5f) / tools::evaluate_polynomial(Q, z - 4.5f);
         result *= exp(-z * z) / z;
      }
      else if(z < 7.5)
      {
         // Maximum Deviation Found:                     1.007e-36
         // Expected Error Term:                         1.007e-36
         // Maximum Relative Change in Control Points:   1.027e-03
         // Max Error found at long double precision =   2.646420e-36
         static const T Y = 0.5574436187744140625f;
         static const T P[] = {    
            0.000293236907400849056269309713064107674L,
            0.00225110719535060642692275221961480162L,
            0.00190984458121502831421717207849429799L,
            0.000747757733460111743833929141001680706L,
            0.000170663175280949889583158597373928096L,
            0.246441188958013822253071608197514058e-4L,
            0.229818000860544644974205957895688106e-5L,
            0.134886977703388748488480980637704864e-6L,
            0.454764611880548962757125070106650958e-8L,
            0.673002744115866600294723141176820155e-10L,
         };
         static const T Q[] = {    
            1L,
            1.12843690320861239631195353379313367L,
            0.569900657061622955362493442186537259L,
            0.169094404206844928112348730277514273L,
            0.0324887449084220415058158657252147063L,
            0.00419252877436825753042680842608219552L,
            0.00036344133176118603523976748563178578L,
            0.204123895931375107397698245752850347e-4L,
            0.674128352521481412232785122943508729e-6L,
            0.997637501418963696542159244436245077e-8L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 6.5f) / tools::evaluate_polynomial(Q, z - 6.5f);
         result *= exp(-z * z) / z;
      }
      else if(z < 11.5)
      {
         // Maximum Deviation Found:                     8.380e-36
         // Expected Error Term:                         8.380e-36
         // Maximum Relative Change in Control Points:   2.632e-06
         // Max Error found at long double precision =   9.849522e-36
         static const T Y = 0.56083202362060546875f;
         static const T P[] = {    
            0.000282420728751494363613829834891390121L,
            0.00175387065018002823433704079355125161L,
            0.0021344978564889819420775336322920375L,
            0.00124151356560137532655039683963075661L,
            0.000423600733566948018555157026862139644L,
            0.914030340865175237133613697319509698e-4L,
            0.126999927156823363353809747017945494e-4L,
            0.110610959842869849776179749369376402e-5L,
            0.55075079477173482096725348704634529e-7L,
            0.119735694018906705225870691331543806e-8L,
         };
         static const T Q[] = {    
            1L,
            1.69889613396167354566098060039549882L,
            1.28824647372749624464956031163282674L,
            0.572297795434934493541628008224078717L,
            0.164157697425571712377043857240773164L,
            0.0315311145224594430281219516531649562L,
            0.00405588922155632380812945849777127458L,
            0.000336929033691445666232029762868642417L,
            0.164033049810404773469413526427932109e-4L,
            0.356615210500531410114914617294694857e-6L,
         };
         result = Y + tools::evaluate_polynomial(P, z / 2 - 4.75f) / tools::evaluate_polynomial(Q, z / 2 - 4.75f);
         result *= exp(-z * z) / z;
      }
      else
      {
         // Maximum Deviation Found:                     1.132e-35
         // Expected Error Term:                         -1.132e-35
         // Maximum Relative Change in Control Points:   4.674e-04
         // Max Error found at long double precision =   1.162590e-35
         static const T Y = 0.5632686614990234375f;
         static const T P[] = {    
            0.000920922048732849448079451574171836943L,
            0.00321439044532288750501700028748922439L,
            -0.250455263029390118657884864261823431L,
            -0.906807635364090342031792404764598142L,
            -8.92233572835991735876688745989985565L,
            -21.7797433494422564811782116907878495L,
            -91.1451915251976354349734589601171659L,
            -144.1279109655993927069052125017673L,
            -313.845076581796338665519022313775589L,
            -273.11378811923343424081101235736475L,
            -271.651566205951067025696102600443452L,
            -60.0530577077238079968843307523245547L,
         };
         static const T Q[] = {    
            1L,
            3.49040448075464744191022350947892036L,
            34.3563592467165971295915749548313227L,
            84.4993232033879023178285731843850461L,
            376.005865281206894120659401340373818L,
            629.95369438888946233003926191755125L,
            1568.35771983533158591604513304269098L,
            1646.02452040831961063640827116581021L,
            2299.96860633240298708910425594484895L,
            1222.73204392037452750381340219906374L,
            799.359797306084372350264298361110448L,
            72.7415265778588087243442792401576737L,
         };
         result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
         result *= exp(-z * z) / z;
      }
   }
   else
   {
      //
      // Any value of z larger than 110 will underflow to zero:
      //
      result = 0;
      invert = !invert;
   }

   if(invert)
   {
      result = 1 - result;
   }

   return result;
} // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<113>& t)

} // namespace detail

template <class T, class Policy>
inline typename tools::promote_args<T>::type erf(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;

   BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name());
   BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name());
   BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name());

   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;

   BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());

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

template <class T, class Policy>
inline typename tools::promote_args<T>::type erfc(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;

   BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name());
   BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name());
   BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name());

   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;

   BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());

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

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

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

} // namespace math
} // namespace boost

#include <boost/math/special_functions/detail/erf_inv.hpp>

#endif // BOOST_MATH_SPECIAL_ERF_HPP