erf.hpp

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

template <class T, class Policy>
T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
{
   BOOST_MATH_STD_USING

   BOOST_MATH_INSTRUMENT_CODE("64-bit precision erf_imp called");

   if(z < 0)
   {
      if(!invert)
         return -erf_imp(-z, invert, pol, t);
      else if(z < -0.5)
         return 2 - erf_imp(-z, invert, pol, t);
      else
         return 1 + erf_imp(-z, false, pol, t);
   }

   T result;

   //
   // Big bunch of selection statements now to pick which
   // implementation to use, try to put most likely options
   // first:
   //
   if(z < 0.5)
   {
      //
      // We're going to calculate erf:
      //
      if(z == 0)
      {
         result = 0;
      }
      else if(z < 1e-10)
      {
         result = z * 1.125 + z * 0.003379167095512573896158903121545171688L;
      }
      else
      {
         // Max Error found at long double precision =   1.623299e-20
         // Maximum Deviation Found:                     4.326e-22
         // Expected Error Term:                         -4.326e-22
         // Maximum Relative Change in Control Points:   1.474e-04
         static const T Y = 1.044948577880859375f;
         static const T P[] = {    
            0.0834305892146531988966L,
            -0.338097283075565413695L,
            -0.0509602734406067204596L,
            -0.00904906346158537794396L,
            -0.000489468651464798669181L,
            -0.200305626366151877759e-4L,
         };
         static const T Q[] = {    
            1L,
            0.455817300515875172439L,
            0.0916537354356241792007L,
            0.0102722652675910031202L,
            0.000650511752687851548735L,
            0.189532519105655496778e-4L,
         };
         result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
      }
   }
   else if((z < 110) || ((z < 110) && invert))  // TODO FIXME!!!
   {
      //
      // We'll be calculating erfc:
      //
      invert = !invert;
      if(z < 1.5)
      {
         // Max Error found at long double precision =   3.239590e-20
         // Maximum Deviation Found:                     2.241e-20
         // Expected Error Term:                         -2.241e-20
         // Maximum Relative Change in Control Points:   5.110e-03
         static const T Y = 0.405935764312744140625f;
         static const T P[] = {    
            -0.0980905922162812031672L,
            0.159989089922969141329L,
            0.222359821619935712378L,
            0.127303921703577362312L,
            0.0384057530342762400273L,
            0.00628431160851156719325L,
            0.000441266654514391746428L,
            0.266689068336295642561e-7L,
         };
         static const T Q[] = {    
            1L,
            2.03237474985469469291L,
            1.78355454954969405222L,
            0.867940326293760578231L,
            0.248025606990021698392L,
            0.0396649631833002269861L,
            0.00279220237309449026796L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
         result *= exp(-z * z) / z;
      }
      else if(z < 2.5)
      {
         // Max Error found at long double precision =   3.686211e-21
         // Maximum Deviation Found:                     1.495e-21
         // Expected Error Term:                         -1.494e-21
         // Maximum Relative Change in Control Points:   1.793e-04
         static const T Y = 0.50672817230224609375f;
         static const T P[] = {    
            -0.024350047620769840217L,
            0.0343522687935671451309L,
            0.0505420824305544949541L,
            0.0257479325917757388209L,
            0.00669349844190354356118L,
            0.00090807914416099524444L,
            0.515917266698050027934e-4L,
         };
         static const T Q[] = {    
            1L,
            1.71657861671930336344L,
            1.26409634824280366218L,
            0.512371437838969015941L,
            0.120902623051120950935L,
            0.0158027197831887485261L,
            0.000897871370778031611439L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
         result *= exp(-z * z) / z;
      }
      else if(z < 4.5)
      {
         // Maximum Deviation Found:                     1.107e-20
         // Expected Error Term:                         -1.106e-20
         // Maximum Relative Change in Control Points:   1.709e-04
         // Max Error found at long double precision =   1.446908e-20
         static const T Y  = 0.5405750274658203125f;
         static const T P[] = {    
            0.0029527671653097284033L,
            0.0141853245895495604051L,
            0.0104959584626432293901L,
            0.00343963795976100077626L,
            0.00059065441194877637899L,
            0.523435380636174008685e-4L,
            0.189896043050331257262e-5L,
         };
         static const T Q[] = {    
            1L,
            1.19352160185285642574L,
            0.603256964363454392857L,
            0.165411142458540585835L,
            0.0259729870946203166468L,
            0.00221657568292893699158L,
            0.804149464190309799804e-4L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 3.5f) / tools::evaluate_polynomial(Q, z - 3.5f);
         result *= exp(-z * z) / z;
      }
      else
      {
         // Max Error found at long double precision =   7.961166e-21
         // Maximum Deviation Found:                     6.677e-21
         // Expected Error Term:                         6.676e-21
         // Maximum Relative Change in Control Points:   2.319e-05
         static const T Y = 0.55825519561767578125f;
         static const T P[] = {    
            0.00593438793008050214106L,
            0.0280666231009089713937L,
            -0.141597835204583050043L,
            -0.978088201154300548842L,
            -5.47351527796012049443L,
            -13.8677304660245326627L,
            -27.1274948720539821722L,
            -29.2545152747009461519L,
            -16.8865774499799676937L,
         };
         static const T Q[] = {    
            1L,
            4.72948911186645394541L,
            23.6750543147695749212L,
            60.0021517335693186785L,
            131.766251645149522868L,
            178.167924971283482513L,
            182.499390505915222699L,
            104.365251479578577989L,
            30.8365511891224291717L,
         };
         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_<64>& t)


template <class T, class Policy>
T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
{
   BOOST_MATH_STD_USING

   BOOST_MATH_INSTRUMENT_CODE("113-bit precision erf_imp called");

   if(z < 0)
   {
      if(!invert)
         return -erf_imp(-z, invert, pol, t);
      else if(z < -0.5)
         return 2 - erf_imp(-z, invert, pol, t);
      else
         return 1 + erf_imp(-z, false, pol, t);
   }

   T result;

   //
   // Big bunch of selection statements now to pick which
   // implementation to use, try to put most likely options
   // first:
   //
   if(z < 0.5)
   {
      //
      // We're going to calculate erf:
      //
      if(z == 0)
      {
         result = 0;
      }
      else if(z < 1e-20)
      {
         result = z * 1.125 + z * 0.003379167095512573896158903121545171688L;
      }
      else
      {
         // Max Error found at long double precision =   2.342380e-35
         // Maximum Deviation Found:                     6.124e-36
         // Expected Error Term:                         -6.124e-36
         // Maximum Relative Change in Control Points:   3.492e-10
         static const T Y = 1.0841522216796875f;
         static const T P[] = {    
            0.0442269454158250738961589031215451778L,
            -0.35549265736002144875335323556961233L,
            -0.0582179564566667896225454670863270393L,
            -0.0112694696904802304229950538453123925L,
            -0.000805730648981801146251825329609079099L,
            -0.566304966591936566229702842075966273e-4L,
            -0.169655010425186987820201021510002265e-5L,
            -0.344448249920445916714548295433198544e-7L,
         };
         static const T Q[] = {    
            1L,
            0.466542092785657604666906909196052522L,
            0.100005087012526447295176964142107611L,
            0.0128341535890117646540050072234142603L,
            0.00107150448466867929159660677016658186L,
            0.586168368028999183607733369248338474e-4L,
            0.196230608502104324965623171516808796e-5L,
            0.313388521582925207734229967907890146e-7L,
         };
         result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
      }
   }
   else if((z < 9) || ((z < 110) && invert))
   {
      //
      // We'll be calculating erfc:
      //
      invert = !invert;
      if(z < 1)
      {
         // Max Error found at long double precision =   3.246278e-35
         // Maximum Deviation Found:                     1.388e-35
         // Expected Error Term:                         1.387e-35
         // Maximum Relative Change in Control Points:   6.127e-05
         static const T Y = 0.371877193450927734375f;
         static const T P[] = {    
            -0.0640320213544647969396032886581290455L,
            0.200769874440155895637857443946706731L,
            0.378447199873537170666487408805779826L,
            0.30521399466465939450398642044975127L,
            0.146890026406815277906781824723458196L,
            0.0464837937749539978247589252732769567L,
            0.00987895759019540115099100165904822903L,
            0.00137507575429025512038051025154301132L,
            0.0001144764551085935580772512359680516L,
            0.436544865032836914773944382339900079e-5L,
         };
         static const T Q[] = {    
            1L,
            2.47651182872457465043733800302427977L,
            2.78706486002517996428836400245547955L,
            1.87295924621659627926365005293130693L,
            0.829375825174365625428280908787261065L,
            0.251334771307848291593780143950311514L,
            0.0522110268876176186719436765734722473L,
            0.00718332151250963182233267040106902368L,
            0.000595279058621482041084986219276392459L,
            0.226988669466501655990637599399326874e-4L,
            0.270666232259029102353426738909226413e-10L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
         result *= exp(-z * z) / z;
      }
      else if(z < 1.5)
      {
         // Max Error found at long double precision =   2.215785e-35
         // Maximum Deviation Found:                     1.539e-35
         // Expected Error Term:                         1.538e-35
         // Maximum Relative Change in Control Points:   6.104e-05
         static const T Y = 0.45658016204833984375f;
         static const T P[] = {    
            -0.0289965858925328393392496555094848345L,
            0.0868181194868601184627743162571779226L,
            0.169373435121178901746317404936356745L,
            0.13350446515949251201104889028133486L,
            0.0617447837290183627136837688446313313L,
            0.0185618495228251406703152962489700468L,
            0.00371949406491883508764162050169531013L,
            0.000485121708792921297742105775823900772L,
            0.376494706741453489892108068231400061e-4L,
            0.133166058052466262415271732172490045e-5L,
         };
         static const T Q[] = {    
            1L,
            2.32970330146503867261275580968135126L,
            2.46325715420422771961250513514928746L,
            1.55307882560757679068505047390857842L,
            0.644274289865972449441174485441409076L,
            0.182609091063258208068606847453955649L,
            0.0354171651271241474946129665801606795L,
            0.00454060370165285246451879969534083997L,
            0.000349871943711566546821198612518656486L,
            0.123749319840299552925421880481085392e-4L,
         };
         result = Y + tools::evaluate_polynomial(P, z - 1.0f) / tools::evaluate_polynomial(Q, z - 1.0f);
         result *= exp(-z * z) / z;
      }
      else if(z < 2.25)
      {
         // Maximum Deviation Found:                     1.418e-35
         // Expected Error Term:                         1.418e-35
         // Maximum Relative Change in Control Points:   1.316e-04
         // Max Error found at long double precision =   1.998462e-35
         static const T Y = 0.50250148773193359375f;
         static const T P[] = {    
            -0.0201233630504573402185161184151016606L,
            0.0331864357574860196516686996302305002L,
            0.0716562720864787193337475444413405461L,
            0.0545835322082103985114927569724880658L,
            0.0236692635189696678976549720784989593L,