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