boost/math/special_functions/erf.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SPECIAL_ERF_HPP
   7 #define BOOST_MATH_SPECIAL_ERF_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/tools/config.hpp>
  15 #include <boost/math/special_functions/gamma.hpp>
  16 #include <boost/math/tools/roots.hpp>
  17 #include <boost/math/policies/error_handling.hpp>
  18
  19 namespace boost{ namespace math{
  20
  21 namespace detail
  22 {
  23
  24 //
  25 // Asymptotic series for large z:
  26 //
  27 template <class T>
  28 struct erf_asympt_series_t
  29 {
  30    erf_asympt_series_t(T z) : xx(2 * -z * z), tk(1)
  31    {
  32       BOOST_MATH_STD_USING
  33       result = -exp(-z * z) / sqrt(boost::math::constants::pi<T>());
  34       result /= z;
  35    }
  36
  37    typedef T result_type;
  38
  39    T operator()()
  40    {
  41       BOOST_MATH_STD_USING
  42       T r = result;
  43       result *= tk / xx;
  44       tk += 2;
  45       if( fabs(r) < fabs(result))
  46          result = 0;
  47       return r;
  48    }
  49 private:
  50    T result;
  51    T xx;
  52    int tk;
  53 };
  54 //
  55 // How large z has to be in order to ensure that the series converges:
  56 //
  57 template <class T>
  58 inline float erf_asymptotic_limit_N(const T&)
  59 {
  60    return (std::numeric_limits<float>::max)();
  61 }
  62 inline float erf_asymptotic_limit_N(const mpl::int_<24>&)
  63 {
  64    return 2.8F;
  65 }
  66 inline float erf_asymptotic_limit_N(const mpl::int_<53>&)
  67 {
  68    return 4.3F;
  69 }
  70 inline float erf_asymptotic_limit_N(const mpl::int_<64>&)
  71 {
  72    return 4.8F;
  73 }
  74 inline float erf_asymptotic_limit_N(const mpl::int_<106>&)
  75 {
  76    return 6.5F;
  77 }
  78 inline float erf_asymptotic_limit_N(const mpl::int_<113>&)
  79 {
  80    return 6.8F;
  81 }
  82
  83 template <class T, class Policy>
  84 inline T erf_asymptotic_limit()
  85 {
  86    typedef typename policies::precision<T, Policy>::type precision_type;
  87    typedef typename mpl::if_<
  88       mpl::less_equal<precision_type, mpl::int_<24> >,
  89       typename mpl::if_<
  90          mpl::less_equal<precision_type, mpl::int_<0> >,
  91          mpl::int_<0>,
  92          mpl::int_<24>
  93       >::type,
  94       typename mpl::if_<
  95          mpl::less_equal<precision_type, mpl::int_<53> >,
  96          mpl::int_<53>,
  97          typename mpl::if_<
  98             mpl::less_equal<precision_type, mpl::int_<64> >,
  99             mpl::int_<64>,
 100             typename mpl::if_<
 101                mpl::less_equal<precision_type, mpl::int_<106> >,
 102                mpl::int_<106>,
 103                typename mpl::if_<
 104                   mpl::less_equal<precision_type, mpl::int_<113> >,
 105                   mpl::int_<113>,
 106                   mpl::int_<0>
 107                >::type
 108             >::type
 109          >::type
 110       >::type
 111    >::type tag_type;
 112    return erf_asymptotic_limit_N(tag_type());
 113 }
 114
 115 template <class T, class Policy, class Tag>
 116 T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
 117 {
 118    BOOST_MATH_STD_USING
 119
 120    BOOST_MATH_INSTRUMENT_CODE("Generic erf_imp called");
 121
 122    if(z < 0)
 123    {
 124       if(!invert)
 125          return -erf_imp(T(-z), invert, pol, t);
 126       else
 127          return 1 + erf_imp(T(-z), false, pol, t);
 128    }
 129
 130    T result;
 131
 132    if(!invert && (z > detail::erf_asymptotic_limit<T, Policy>()))
 133    {
 134       detail::erf_asympt_series_t<T> s(z);
 135       boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 136       result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, 1);
 137       policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
 138    }
 139    else
 140    {
 141       T x = z * z;
 142       if(x < 0.6)
 143       {
 144          // Compute P:
 145          result = z * exp(-x);
 146          result /= sqrt(boost::math::constants::pi<T>());
 147          if(result != 0)
 148             result *= 2 * detail::lower_gamma_series(T(0.5f), x, pol);
 149       }
 150       else if(x < 1.1f)
 151       {
 152          // Compute Q:
 153          invert = !invert;
 154          result = tgamma_small_upper_part(T(0.5f), x, pol);
 155          result /= sqrt(boost::math::constants::pi<T>());
 156       }
 157       else
 158       {
 159          // Compute Q:
 160          invert = !invert;
 161          result = z * exp(-x);
 162          result /= sqrt(boost::math::constants::pi<T>());
 163          result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon<T, Policy>());
 164       }
 165    }
 166    if(invert)
 167       result = 1 - result;
 168    return result;
 169 }
 170
 171 template <class T, class Policy>
 172 T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
 173 {
 174    BOOST_MATH_STD_USING
 175
 176    BOOST_MATH_INSTRUMENT_CODE("53-bit precision erf_imp called");
 177
 178    if(z < 0)
 179    {
 180       if(!invert)
 181          return -erf_imp(-z, invert, pol, t);
 182       else if(z < -0.5)
 183          return 2 - erf_imp(-z, invert, pol, t);
 184       else
 185          return 1 + erf_imp(-z, false, pol, t);
 186    }
 187
 188    T result;
 189
 190    //
 191    // Big bunch of selection statements now to pick
 192    // which implementation to use,
 193    // try to put most likely options first:
 194    //
 195    if(z < 0.5)
 196    {
 197       //
 198       // We're going to calculate erf:
 199       //
 200       if(z < 1e-10)
 201       {
 202          if(z == 0)
 203          {
 204             result = T(0);
 205          }
 206          else
 207          {
 208             result = static_cast<T>(z * 1.125f + z * 0.003379167095512573896158903121545171688L);
 209          }
 210       }
 211       else
 212       {
 213          // Maximum Deviation Found:                     1.561e-17
 214          // Expected Error Term:                         1.561e-17
 215          // Maximum Relative Change in Control Points:   1.155e-04
 216          // Max Error found at double precision =        2.961182e-17
 217
 218          static const T Y = 1.044948577880859375f;
 219          static const T P[] = {
 220             0.0834305892146531832907L,
 221             -0.338165134459360935041L,
 222             -0.0509990735146777432841L,
 223             -0.00772758345802133288487L,
 224             -0.000322780120964605683831L,
 225          };
 226          static const T Q[] = {
 227             1L,
 228             0.455004033050794024546L,
 229             0.0875222600142252549554L,
 230             0.00858571925074406212772L,
 231             0.000370900071787748000569L,
 232          };
 233          T zz = z * z;
 234          result = z * (Y + tools::evaluate_polynomial(P, zz) / tools::evaluate_polynomial(Q, zz));
 235       }
 236    }
 237    else if(invert ? (z < 28) : (z < 5.8f))
 238    {
 239       //
 240       // We'll be calculating erfc:
 241       //
 242       invert = !invert;
 243       if(z < 1.5f)
 244       {
 245          // Maximum Deviation Found:                     3.702e-17
 246          // Expected Error Term:                         3.702e-17
 247          // Maximum Relative Change in Control Points:   2.845e-04
 248          // Max Error found at double precision =        4.841816e-17
 249          static const T Y = 0.405935764312744140625f;
 250          static const T P[] = {
 251             -0.098090592216281240205L,
 252             0.178114665841120341155L,
 253             0.191003695796775433986L,
 254             0.0888900368967884466578L,
 255             0.0195049001251218801359L,
 256             0.00180424538297014223957L,
 257          };
 258          static const T Q[] = {
 259             1L,
 260             1.84759070983002217845L,
 261             1.42628004845511324508L,
 262             0.578052804889902404909L,
 263             0.12385097467900864233L,
 264             0.0113385233577001411017L,
 265             0.337511472483094676155e-5L,
 266          };
 267          result = Y + tools::evaluate_polynomial(P, z - 0.5) / tools::evaluate_polynomial(Q, z - 0.5);
 268          result *= exp(-z * z) / z;
 269       }
 270       else if(z < 2.5f)
 271       {
 272          // Max Error found at double precision =        6.599585e-18
 273          // Maximum Deviation Found:                     3.909e-18
 274          // Expected Error Term:                         3.909e-18
 275          // Maximum Relative Change in Control Points:   9.886e-05
 276          static const T Y = 0.50672817230224609375f;
 277          static const T P[] = {
 278             -0.0243500476207698441272L,
 279             0.0386540375035707201728L,
 280             0.04394818964209516296L,
 281             0.0175679436311802092299L,
 282             0.00323962406290842133584L,
 283             0.000235839115596880717416L,
 284          };
 285          static const T Q[] = {
 286             1L,
 287             1.53991494948552447182L,
 288             0.982403709157920235114L,
 289             0.325732924782444448493L,
 290             0.0563921837420478160373L,
 291             0.00410369723978904575884L,
 292          };
 293          result = Y + tools::evaluate_polynomial(P, z - 1.5) / tools::evaluate_polynomial(Q, z - 1.5);
 294          result *= exp(-z * z) / z;
 295       }
 296       else if(z < 4.5f)
 297       {
 298          // Maximum Deviation Found:                     1.512e-17
 299          // Expected Error Term:                         1.512e-17
 300          // Maximum Relative Change in Control Points:   2.222e-04
 301          // Max Error found at double precision =        2.062515e-17
 302          static const T Y = 0.5405750274658203125f;
 303          static const T P[] = {
 304             0.00295276716530971662634L,
 305             0.0137384425896355332126L,
 306             0.00840807615555585383007L,
 307             0.00212825620914618649141L,
 308             0.000250269961544794627958L,
 309             0.113212406648847561139e-4L,
 310          };
 311          static const T Q[] = {
 312             1L,
 313             1.04217814166938418171L,
 314             0.442597659481563127003L,
 315             0.0958492726301061423444L,
 316             0.0105982906484876531489L,
 317             0.000479411269521714493907L,
 318          };
 319          result = Y + tools::evaluate_polynomial(P, z - 3.5) / tools::evaluate_polynomial(Q, z - 3.5);
 320          result *= exp(-z * z) / z;
 321       }
 322       else
 323       {
 324          // Max Error found at double precision =        2.997958e-17
 325          // Maximum Deviation Found:                     2.860e-17
 326          // Expected Error Term:                         2.859e-17
 327          // Maximum Relative Change in Control Points:   1.357e-05
 328          static const T Y = 0.5579090118408203125f;
 329          static const T P[] = {
 330             0.00628057170626964891937L,
 331             0.0175389834052493308818L,
 332             -0.212652252872804219852L,
 333             -0.687717681153649930619L,
 334             -2.5518551727311523996L,
 335             -3.22729451764143718517L,
 336             -2.8175401114513378771L,
 337          };
 338          static const T Q[] = {
 339             1L,
 340             2.79257750980575282228L,
 341             11.0567237927800161565L,
 342             15.930646027911794143L,
 343             22.9367376522880577224L,
 344             13.5064170191802889145L,
 345             5.48409182238641741584L,
 346          };
 347          result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
 348          result *= exp(-z * z) / z;
 349       }
 350    }
 351    else
 352    {
 353       //
 354       // Any value of z larger than 28 will underflow to zero:
 355       //
 356       result = 0;
 357       invert = !invert;
 358    }
 359
 360    if(invert)
 361    {
 362       result = 1 - result;
 363    }
 364
 365    return result;
 366 } // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<53>& t)
 367
 368
 369 template <class T, class Policy>
 370 T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
 371 {
 372    BOOST_MATH_STD_USING
 373
 374    BOOST_MATH_INSTRUMENT_CODE("64-bit precision erf_imp called");
 375
 376    if(z < 0)
 377    {
 378       if(!invert)
 379          return -erf_imp(-z, invert, pol, t);
 380       else if(z < -0.5)
 381          return 2 - erf_imp(-z, invert, pol, t);
 382       else
 383          return 1 + erf_imp(-z, false, pol, t);
 384    }
 385
 386    T result;
 387
 388    //
 389    // Big bunch of selection statements now to pick which
 390    // implementation to use, try to put most likely options
 391    // first:
 392    //
 393    if(z < 0.5)
 394    {
 395       //
 396       // We're going to calculate erf:
 397       //
 398       if(z == 0)
 399       {
 400          result = 0;
 401       }
 402       else if(z < 1e-10)
 403       {
 404          result = z * 1.125 + z * 0.003379167095512573896158903121545171688L;
 405       }
 406       else
 407       {
 408          // Max Error found at long double precision =   1.623299e-20
 409          // Maximum Deviation Found:                     4.326e-22
 410          // Expected Error Term:                         -4.326e-22
 411          // Maximum Relative Change in Control Points:   1.474e-04
 412          static const T Y = 1.044948577880859375f;
 413          static const T P[] = {
 414             0.0834305892146531988966L,
 415             -0.338097283075565413695L,
 416             -0.0509602734406067204596L,
 417             -0.00904906346158537794396L,
 418             -0.000489468651464798669181L,
 419             -0.200305626366151877759e-4L,
 420          };
 421          static const T Q[] = {
 422             1L,
 423             0.455817300515875172439L,
 424             0.0916537354356241792007L,
 425             0.0102722652675910031202L,
 426             0.000650511752687851548735L,
 427             0.189532519105655496778e-4L,
 428          };
 429          result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
 430       }
 431    }
 432    else if(invert ? (z < 110) : (z < 6.4f))
 433    {
 434       //
 435       // We'll be calculating erfc:
 436       //
 437       invert = !invert;
 438       if(z < 1.5)
 439       {
 440          // Max Error found at long double precision =   3.239590e-20
 441          // Maximum Deviation Found:                     2.241e-20
 442          // Expected Error Term:                         -2.241e-20
 443          // Maximum Relative Change in Control Points:   5.110e-03
 444          static const T Y = 0.405935764312744140625f;
 445          static const T P[] = {
 446             -0.0980905922162812031672L,
 447             0.159989089922969141329L,
 448             0.222359821619935712378L,
 449             0.127303921703577362312L,
 450             0.0384057530342762400273L,
 451             0.00628431160851156719325L,
 452             0.000441266654514391746428L,
 453             0.266689068336295642561e-7L,
 454          };
 455          static const T Q[] = {
 456             1L,
 457             2.03237474985469469291L,
 458             1.78355454954969405222L,
 459             0.867940326293760578231L,
 460             0.248025606990021698392L,
 461             0.0396649631833002269861L,
 462             0.00279220237309449026796L,
 463          };
 464          result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
 465          result *= exp(-z * z) / z;
 466       }
 467       else if(z < 2.5)
 468       {
 469          // Max Error found at long double precision =   3.686211e-21
 470          // Maximum Deviation Found:                     1.495e-21
 471          // Expected Error Term:                         -1.494e-21
 472          // Maximum Relative Change in Control Points:   1.793e-04
 473          static const T Y = 0.50672817230224609375f;
 474          static const T P[] = {
 475             -0.024350047620769840217L,
 476             0.0343522687935671451309L,
 477             0.0505420824305544949541L,
 478             0.0257479325917757388209L,
 479             0.00669349844190354356118L,
 480             0.00090807914416099524444L,
 481             0.515917266698050027934e-4L,
 482          };
 483          static const T Q[] = {
 484             1L,
 485             1.71657861671930336344L,
 486             1.26409634824280366218L,
 487             0.512371437838969015941L,
 488             0.120902623051120950935L,
 489             0.0158027197831887485261L,
 490             0.000897871370778031611439L,
 491          };
 492          result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
 493          result *= exp(-z * z) / z;
 494       }
 495       else if(z < 4.5)
 496       {
 497          // Maximum Deviation Found:                     1.107e-20
 498          // Expected Error Term:                         -1.106e-20
 499          // Maximum Relative Change in Control Points:   1.709e-04
 500          // Max Error found at long double precision =   1.446908e-20
 501          static const T Y  = 0.5405750274658203125f;
 502          static const T P[] = {
 503             0.0029527671653097284033L,
 504             0.0141853245895495604051L,
 505             0.0104959584626432293901L,
 506             0.00343963795976100077626L,
 507             0.00059065441194877637899L,
 508             0.523435380636174008685e-4L,
 509             0.189896043050331257262e-5L,
 510          };
 511          static const T Q[] = {
 512             1L,
 513             1.19352160185285642574L,
 514             0.603256964363454392857L,
 515             0.165411142458540585835L,
 516             0.0259729870946203166468L,
 517             0.00221657568292893699158L,
 518             0.804149464190309799804e-4L,
 519          };
 520          result = Y + tools::evaluate_polynomial(P, z - 3.5f) / tools::evaluate_polynomial(Q, z - 3.5f);
 521          result *= exp(-z * z) / z;
 522       }
 523       else
 524       {
 525          // Max Error found at long double precision =   7.961166e-21
 526          // Maximum Deviation Found:                     6.677e-21
 527          // Expected Error Term:                         6.676e-21
 528          // Maximum Relative Change in Control Points:   2.319e-05
 529          static const T Y = 0.55825519561767578125f;
 530          static const T P[] = {
 531             0.00593438793008050214106L,
 532             0.0280666231009089713937L,
 533             -0.141597835204583050043L,
 534             -0.978088201154300548842L,
 535             -5.47351527796012049443L,
 536             -13.8677304660245326627L,
 537             -27.1274948720539821722L,
 538             -29.2545152747009461519L,
 539             -16.8865774499799676937L,
 540          };
 541          static const T Q[] = {
 542             1L,
 543             4.72948911186645394541L,
 544             23.6750543147695749212L,
 545             60.0021517335693186785L,
 546             131.766251645149522868L,
 547             178.167924971283482513L,
 548             182.499390505915222699L,
 549             104.365251479578577989L,
 550             30.8365511891224291717L,
 551          };
 552          result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
 553          result *= exp(-z * z) / z;
 554       }
 555    }
 556    else
 557    {
 558       //
 559       // Any value of z larger than 110 will underflow to zero:
 560       //
 561       result = 0;
 562       invert = !invert;
 563    }
 564
 565    if(invert)
 566    {
 567       result = 1 - result;
 568    }
 569
 570    return result;
 571 } // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<64>& t)
 572
 573
 574 template <class T, class Policy>
 575 T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
 576 {
 577    BOOST_MATH_STD_USING
 578
 579    BOOST_MATH_INSTRUMENT_CODE("113-bit precision erf_imp called");
 580
 581    if(z < 0)
 582    {
 583       if(!invert)
 584          return -erf_imp(-z, invert, pol, t);
 585       else if(z < -0.5)
 586          return 2 - erf_imp(-z, invert, pol, t);
 587       else
 588          return 1 + erf_imp(-z, false, pol, t);
 589    }
 590
 591    T result;
 592
 593    //
 594    // Big bunch of selection statements now to pick which
 595    // implementation to use, try to put most likely options
 596    // first:
 597    //
 598    if(z < 0.5)
 599    {
 600       //
 601       // We're going to calculate erf:
 602       //
 603       if(z == 0)
 604       {
 605          result = 0;
 606       }
 607       else if(z < 1e-20)
 608       {
 609          result = z * 1.125 + z * 0.003379167095512573896158903121545171688L;
 610       }
 611       else
 612       {
 613          // Max Error found at long double precision =   2.342380e-35
 614          // Maximum Deviation Found:                     6.124e-36
 615          // Expected Error Term:                         -6.124e-36
 616          // Maximum Relative Change in Control Points:   3.492e-10
 617          static const T Y = 1.0841522216796875f;
 618          static const T P[] = {
 619             0.0442269454158250738961589031215451778L,
 620             -0.35549265736002144875335323556961233L,
 621             -0.0582179564566667896225454670863270393L,
 622             -0.0112694696904802304229950538453123925L,
 623             -0.000805730648981801146251825329609079099L,
 624             -0.566304966591936566229702842075966273e-4L,
 625             -0.169655010425186987820201021510002265e-5L,
 626             -0.344448249920445916714548295433198544e-7L,
 627          };
 628          static const T Q[] = {
 629             1L,
 630             0.466542092785657604666906909196052522L,
 631             0.100005087012526447295176964142107611L,
 632             0.0128341535890117646540050072234142603L,
 633             0.00107150448466867929159660677016658186L,
 634             0.586168368028999183607733369248338474e-4L,
 635             0.196230608502104324965623171516808796e-5L,
 636             0.313388521582925207734229967907890146e-7L,
 637          };
 638          result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
 639       }
 640    }
 641    else if(invert ? (z < 110) : (z < 8.65f))
 642    {
 643       //
 644       // We'll be calculating erfc:
 645       //
 646       invert = !invert;
 647       if(z < 1)
 648       {
 649          // Max Error found at long double precision =   3.246278e-35
 650          // Maximum Deviation Found:                     1.388e-35
 651          // Expected Error Term:                         1.387e-35
 652          // Maximum Relative Change in Control Points:   6.127e-05
 653          static const T Y = 0.371877193450927734375f;
 654          static const T P[] = {
 655             -0.0640320213544647969396032886581290455L,
 656             0.200769874440155895637857443946706731L,
 657             0.378447199873537170666487408805779826L,
 658             0.30521399466465939450398642044975127L,
 659             0.146890026406815277906781824723458196L,
 660             0.0464837937749539978247589252732769567L,
 661             0.00987895759019540115099100165904822903L,
 662             0.00137507575429025512038051025154301132L,
 663             0.0001144764551085935580772512359680516L,
 664             0.436544865032836914773944382339900079e-5L,
 665          };
 666          static const T Q[] = {
 667             1L,
 668             2.47651182872457465043733800302427977L,
 669             2.78706486002517996428836400245547955L,
 670             1.87295924621659627926365005293130693L,
 671             0.829375825174365625428280908787261065L,
 672             0.251334771307848291593780143950311514L,
 673             0.0522110268876176186719436765734722473L,
 674             0.00718332151250963182233267040106902368L,
 675             0.000595279058621482041084986219276392459L,
 676             0.226988669466501655990637599399326874e-4L,
 677             0.270666232259029102353426738909226413e-10L,
 678          };
 679          result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
 680          result *= exp(-z * z) / z;
 681       }
 682       else if(z < 1.5)
 683       {
 684          // Max Error found at long double precision =   2.215785e-35
 685          // Maximum Deviation Found:                     1.539e-35
 686          // Expected Error Term:                         1.538e-35
 687          // Maximum Relative Change in Control Points:   6.104e-05
 688          static const T Y = 0.45658016204833984375f;
 689          static const T P[] = {
 690             -0.0289965858925328393392496555094848345L,
 691             0.0868181194868601184627743162571779226L,
 692             0.169373435121178901746317404936356745L,
 693             0.13350446515949251201104889028133486L,
 694             0.0617447837290183627136837688446313313L,
 695             0.0185618495228251406703152962489700468L,
 696             0.00371949406491883508764162050169531013L,
 697             0.000485121708792921297742105775823900772L,
 698             0.376494706741453489892108068231400061e-4L,
 699             0.133166058052466262415271732172490045e-5L,
 700          };
 701          static const T Q[] = {
 702             1L,
 703             2.32970330146503867261275580968135126L,
 704             2.46325715420422771961250513514928746L,
 705             1.55307882560757679068505047390857842L,
 706             0.644274289865972449441174485441409076L,
 707             0.182609091063258208068606847453955649L,
 708             0.0354171651271241474946129665801606795L,
 709             0.00454060370165285246451879969534083997L,
 710             0.000349871943711566546821198612518656486L,
 711             0.123749319840299552925421880481085392e-4L,
 712          };
 713          result = Y + tools::evaluate_polynomial(P, z - 1.0f) / tools::evaluate_polynomial(Q, z - 1.0f);
 714          result *= exp(-z * z) / z;
 715       }
 716       else if(z < 2.25)
 717       {
 718          // Maximum Deviation Found:                     1.418e-35
 719          // Expected Error Term:                         1.418e-35
 720          // Maximum Relative Change in Control Points:   1.316e-04
 721          // Max Error found at long double precision =   1.998462e-35
 722          static const T Y = 0.50250148773193359375f;
 723          static const T P[] = {
 724             -0.0201233630504573402185161184151016606L,
 725             0.0331864357574860196516686996302305002L,
 726             0.0716562720864787193337475444413405461L,
 727             0.0545835322082103985114927569724880658L,
 728             0.0236692635189696678976549720784989593L,
 729             0.00656970902163248872837262539337601845L,
 730             0.00120282643299089441390490459256235021L,
 731             0.000142123229065182650020762792081622986L,
 732             0.991531438367015135346716277792989347e-5L,
 733             0.312857043762117596999398067153076051e-6L,
 734          };
 735          static const T Q[] = {
 736             1L,
 737             2.13506082409097783827103424943508554L,
 738             2.06399257267556230937723190496806215L,
 739             1.18678481279932541314830499880691109L,
 740             0.447733186643051752513538142316799562L,
 741             0.11505680005657879437196953047542148L,
 742             0.020163993632192726170219663831914034L,
 743             0.00232708971840141388847728782209730585L,
 744             0.000160733201627963528519726484608224112L,
 745             0.507158721790721802724402992033269266e-5L,
 746             0.18647774409821470950544212696270639e-12L,
 747          };
 748          result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
 749          result *= exp(-z * z) / z;
 750       }
 751       else if (z < 3)
 752       {
 753          // Maximum Deviation Found:                     3.575e-36
 754          // Expected Error Term:                         3.575e-36
 755          // Maximum Relative Change in Control Points:   7.103e-05
 756          // Max Error found at long double precision =   5.794737e-36
 757          static const T Y = 0.52896785736083984375f;
 758          static const T P[] = {
 759             -0.00902152521745813634562524098263360074L,
 760             0.0145207142776691539346923710537580927L,
 761             0.0301681239582193983824211995978678571L,
 762             0.0215548540823305814379020678660434461L,
 763             0.00864683476267958365678294164340749949L,
 764             0.00219693096885585491739823283511049902L,
 765             0.000364961639163319762492184502159894371L,
 766             0.388174251026723752769264051548703059e-4L,
 767             0.241918026931789436000532513553594321e-5L,
 768             0.676586625472423508158937481943649258e-7L,
 769          };
 770          static const T Q[] = {
 771             1L,
 772             1.93669171363907292305550231764920001L,
 773             1.69468476144051356810672506101377494L,
 774             0.880023580986436640372794392579985511L,
 775             0.299099106711315090710836273697708402L,
 776             0.0690593962363545715997445583603382337L,
 777             0.0108427016361318921960863149875360222L,
 778             0.00111747247208044534520499324234317695L,
 779             0.686843205749767250666787987163701209e-4L,
 780             0.192093541425429248675532015101904262e-5L,
 781          };
 782          result = Y + tools::evaluate_polynomial(P, z - 2.25f) / tools::evaluate_polynomial(Q, z - 2.25f);
 783          result *= exp(-z * z) / z;
 784       }
 785       else if(z < 3.5)
 786       {
 787          // Maximum Deviation Found:                     8.126e-37
 788          // Expected Error Term:                         -8.126e-37
 789          // Maximum Relative Change in Control Points:   1.363e-04
 790          // Max Error found at long double precision =   1.747062e-36
 791          static const T Y = 0.54037380218505859375f;
 792          static const T P[] = {
 793             -0.0033703486408887424921155540591370375L,
 794             0.0104948043110005245215286678898115811L,
 795             0.0148530118504000311502310457390417795L,
 796             0.00816693029245443090102738825536188916L,
 797             0.00249716579989140882491939681805594585L,
 798             0.0004655591010047353023978045800916647L,
 799             0.531129557920045295895085236636025323e-4L,
 800             0.343526765122727069515775194111741049e-5L,
 801             0.971120407556888763695313774578711839e-7L,
 802          };
 803          static const T Q[] = {
 804             1L,
 805             1.59911256167540354915906501335919317L,
 806             1.136006830764025173864831382946934L,
 807             0.468565867990030871678574840738423023L,
 808             0.122821824954470343413956476900662236L,
 809             0.0209670914950115943338996513330141633L,
 810             0.00227845718243186165620199012883547257L,
 811             0.000144243326443913171313947613547085553L,
 812             0.407763415954267700941230249989140046e-5L,
 813          };
 814          result = Y + tools::evaluate_polynomial(P, z - 3.0f) / tools::evaluate_polynomial(Q, z - 3.0f);
 815          result *= exp(-z * z) / z;
 816       }
 817       else if(z < 5.5)
 818       {
 819          // Maximum Deviation Found:                     5.804e-36
 820          // Expected Error Term:                         -5.803e-36
 821          // Maximum Relative Change in Control Points:   2.475e-05
 822          // Max Error found at long double precision =   1.349545e-35
 823          static const T Y = 0.55000019073486328125f;
 824          static const T P[] = {
 825             0.00118142849742309772151454518093813615L,
 826             0.0072201822885703318172366893469382745L,
 827             0.0078782276276860110721875733778481505L,
 828             0.00418229166204362376187593976656261146L,
 829             0.00134198400587769200074194304298642705L,
 830             0.000283210387078004063264777611497435572L,
 831             0.405687064094911866569295610914844928e-4L,
 832             0.39348283801568113807887364414008292e-5L,
 833             0.248798540917787001526976889284624449e-6L,
 834             0.929502490223452372919607105387474751e-8L,
 835             0.156161469668275442569286723236274457e-9L,
 836          };
 837          static const T Q[] = {
 838             1L,
 839             1.52955245103668419479878456656709381L,
 840             1.06263944820093830054635017117417064L,
 841             0.441684612681607364321013134378316463L,
 842             0.121665258426166960049773715928906382L,
 843             0.0232134512374747691424978642874321434L,
 844             0.00310778180686296328582860464875562636L,
 845             0.000288361770756174705123674838640161693L,
 846             0.177529187194133944622193191942300132e-4L,
 847             0.655068544833064069223029299070876623e-6L,
 848             0.11005507545746069573608988651927452e-7L,
 849          };
 850          result = Y + tools::evaluate_polynomial(P, z - 4.5f) / tools::evaluate_polynomial(Q, z - 4.5f);
 851          result *= exp(-z * z) / z;
 852       }
 853       else if(z < 7.5)
 854       {
 855          // Maximum Deviation Found:                     1.007e-36
 856          // Expected Error Term:                         1.007e-36
 857          // Maximum Relative Change in Control Points:   1.027e-03
 858          // Max Error found at long double precision =   2.646420e-36
 859          static const T Y = 0.5574436187744140625f;
 860          static const T P[] = {
 861             0.000293236907400849056269309713064107674L,
 862             0.00225110719535060642692275221961480162L,
 863             0.00190984458121502831421717207849429799L,
 864             0.000747757733460111743833929141001680706L,
 865             0.000170663175280949889583158597373928096L,
 866             0.246441188958013822253071608197514058e-4L,
 867             0.229818000860544644974205957895688106e-5L,
 868             0.134886977703388748488480980637704864e-6L,
 869             0.454764611880548962757125070106650958e-8L,
 870             0.673002744115866600294723141176820155e-10L,
 871          };
 872          static const T Q[] = {
 873             1L,
 874             1.12843690320861239631195353379313367L,
 875             0.569900657061622955362493442186537259L,
 876             0.169094404206844928112348730277514273L,
 877             0.0324887449084220415058158657252147063L,
 878             0.00419252877436825753042680842608219552L,
 879             0.00036344133176118603523976748563178578L,
 880             0.204123895931375107397698245752850347e-4L,
 881             0.674128352521481412232785122943508729e-6L,
 882             0.997637501418963696542159244436245077e-8L,
 883          };
 884          result = Y + tools::evaluate_polynomial(P, z - 6.5f) / tools::evaluate_polynomial(Q, z - 6.5f);
 885          result *= exp(-z * z) / z;
 886       }
 887       else if(z < 11.5)
 888       {
 889          // Maximum Deviation Found:                     8.380e-36
 890          // Expected Error Term:                         8.380e-36
 891          // Maximum Relative Change in Control Points:   2.632e-06
 892          // Max Error found at long double precision =   9.849522e-36
 893          static const T Y = 0.56083202362060546875f;
 894          static const T P[] = {
 895             0.000282420728751494363613829834891390121L,
 896             0.00175387065018002823433704079355125161L,
 897             0.0021344978564889819420775336322920375L,
 898             0.00124151356560137532655039683963075661L,
 899             0.000423600733566948018555157026862139644L,
 900             0.914030340865175237133613697319509698e-4L,
 901             0.126999927156823363353809747017945494e-4L,
 902             0.110610959842869849776179749369376402e-5L,
 903             0.55075079477173482096725348704634529e-7L,
 904             0.119735694018906705225870691331543806e-8L,
 905          };
 906          static const T Q[] = {
 907             1L,
 908             1.69889613396167354566098060039549882L,
 909             1.28824647372749624464956031163282674L,
 910             0.572297795434934493541628008224078717L,
 911             0.164157697425571712377043857240773164L,
 912             0.0315311145224594430281219516531649562L,
 913             0.00405588922155632380812945849777127458L,
 914             0.000336929033691445666232029762868642417L,
 915             0.164033049810404773469413526427932109e-4L,
 916             0.356615210500531410114914617294694857e-6L,
 917          };
 918          result = Y + tools::evaluate_polynomial(P, z / 2 - 4.75f) / tools::evaluate_polynomial(Q, z / 2 - 4.75f);
 919          result *= exp(-z * z) / z;
 920       }
 921       else
 922       {
 923          // Maximum Deviation Found:                     1.132e-35
 924          // Expected Error Term:                         -1.132e-35
 925          // Maximum Relative Change in Control Points:   4.674e-04
 926          // Max Error found at long double precision =   1.162590e-35
 927          static const T Y = 0.5632686614990234375f;
 928          static const T P[] = {
 929             0.000920922048732849448079451574171836943L,
 930             0.00321439044532288750501700028748922439L,
 931             -0.250455263029390118657884864261823431L,
 932             -0.906807635364090342031792404764598142L,
 933             -8.92233572835991735876688745989985565L,
 934             -21.7797433494422564811782116907878495L,
 935             -91.1451915251976354349734589601171659L,
 936             -144.1279109655993927069052125017673L,
 937             -313.845076581796338665519022313775589L,
 938             -273.11378811923343424081101235736475L,
 939             -271.651566205951067025696102600443452L,
 940             -60.0530577077238079968843307523245547L,
 941          };
 942          static const T Q[] = {
 943             1L,
 944             3.49040448075464744191022350947892036L,
 945             34.3563592467165971295915749548313227L,
 946             84.4993232033879023178285731843850461L,
 947             376.005865281206894120659401340373818L,
 948             629.95369438888946233003926191755125L,
 949             1568.35771983533158591604513304269098L,
 950             1646.02452040831961063640827116581021L,
 951             2299.96860633240298708910425594484895L,
 952             1222.73204392037452750381340219906374L,
 953             799.359797306084372350264298361110448L,
 954             72.7415265778588087243442792401576737L,
 955          };
 956          result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
 957          result *= exp(-z * z) / z;
 958       }
 959    }
 960    else
 961    {
 962       //
 963       // Any value of z larger than 110 will underflow to zero:
 964       //
 965       result = 0;
 966       invert = !invert;
 967    }
 968
 969    if(invert)
 970    {
 971       result = 1 - result;
 972    }
 973
 974    return result;
 975 } // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<113>& t)
 976
 977 } // namespace detail
 978
 979 template <class T, class Policy>
 980 inline typename tools::promote_args<T>::type erf(T z, const Policy& /* pol */)
 981 {
 982    typedef typename tools::promote_args<T>::type result_type;
 983    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 984    typedef typename policies::precision<result_type, Policy>::type precision_type;
 985    typedef typename policies::normalise<
 986       Policy,
 987       policies::promote_float<false>,
 988       policies::promote_double<false>,
 989       policies::discrete_quantile<>,
 990       policies::assert_undefined<> >::type forwarding_policy;
 991
 992    BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name());
 993    BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name());
 994    BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name());
 995
 996    typedef typename mpl::if_<
 997       mpl::less_equal<precision_type, mpl::int_<0> >,
 998       mpl::int_<0>,
 999       typename mpl::if_<
1000          mpl::less_equal<precision_type, mpl::int_<53> >,
1001          mpl::int_<53>,  // double
1002          typename mpl::if_<
1003             mpl::less_equal<precision_type, mpl::int_<64> >,
1004             mpl::int_<64>, // 80-bit long double
1005             typename mpl::if_<
1006                mpl::less_equal<precision_type, mpl::int_<113> >,
1007                mpl::int_<113>, // 128-bit long double
1008                mpl::int_<0> // too many bits, use generic version.
1009             >::type
1010          >::type
1011       >::type
1012    >::type tag_type;
1013
1014    BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());
1015
1016    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp(
1017       static_cast<value_type>(z),
1018       false,
1019       forwarding_policy(),
1020       tag_type()), "boost::math::erf<%1%>(%1%, %1%)");
1021 }
1022
1023 template <class T, class Policy>
1024 inline typename tools::promote_args<T>::type erfc(T z, const Policy& /* pol */)
1025 {
1026    typedef typename tools::promote_args<T>::type result_type;
1027    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1028    typedef typename policies::precision<result_type, Policy>::type precision_type;
1029    typedef typename policies::normalise<
1030       Policy,
1031       policies::promote_float<false>,
1032       policies::promote_double<false>,
1033       policies::discrete_quantile<>,
1034       policies::assert_undefined<> >::type forwarding_policy;
1035
1036    BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name());
1037    BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name());
1038    BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name());
1039
1040    typedef typename mpl::if_<
1041       mpl::less_equal<precision_type, mpl::int_<0> >,
1042       mpl::int_<0>,
1043       typename mpl::if_<
1044          mpl::less_equal<precision_type, mpl::int_<53> >,
1045          mpl::int_<53>,  // double
1046          typename mpl::if_<
1047             mpl::less_equal<precision_type, mpl::int_<64> >,
1048             mpl::int_<64>, // 80-bit long double
1049             typename mpl::if_<
1050                mpl::less_equal<precision_type, mpl::int_<113> >,
1051                mpl::int_<113>, // 128-bit long double
1052                mpl::int_<0> // too many bits, use generic version.
1053             >::type
1054          >::type
1055       >::type
1056    >::type tag_type;
1057
1058    BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());
1059
1060    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp(
1061       static_cast<value_type>(z),
1062       true,
1063       forwarding_policy(),
1064       tag_type()), "boost::math::erfc<%1%>(%1%, %1%)");
1065 }
1066
1067 template <class T>
1068 inline typename tools::promote_args<T>::type erf(T z)
1069 {
1070    return boost::math::erf(z, policies::policy<>());
1071 }
1072
1073 template <class T>
1074 inline typename tools::promote_args<T>::type erfc(T z)
1075 {
1076    return boost::math::erfc(z, policies::policy<>());
1077 }
1078
1079 } // namespace math
1080 } // namespace boost
1081
1082 #include <boost/math/special_functions/detail/erf_inv.hpp>
1083
1084 #endif // BOOST_MATH_SPECIAL_ERF_HPP
1085
1086
1087
1088